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Chapter  1 


INTRODUCTION  AND  LITERATURE  REVIEW 

During  steady  operation  of  a  liquid  propellant  rocket  engine  the 
injected  propellants  are  converted  by  various  physical  and  chemical 
processes  into  hot  burned  gases  which  are  subsequently  accelerated  to 
supersonic  velocity  by  passing  through  a  converging-diverging  nozzle.  The 
operation  of  such  an  engine,  however,  is  seldom  perfectly  smooth.  Instead 
the  quantities  which  describe  the  conditions  inside  the  combustor  (i.e. 
pressure,  density,  temperature,  etc.)  are  time-dependent  and  oscillatory. 
Such  oscillations  can  be  of  either  a  destructive  or  nondestructive  nature. 
Nondestructive  unsteadiness  is  characterized  by  random  fluctuations  in  the 
flow  properties  and  includes  the  phenomena  of  turbulence  and  combustion 
noise.  Unsteady  operation  of  a  destructive  nature,  on  the  other  hand,  is 
characterized  by  organized  oscillations  in  which  there  is  a  definite 
correlation  between  the  fluctuations  at  two  different  locations  in  the 
combustor.  Such  oscillations  have  a  definite  frequency  and  result  in 
additional  thermal  and  mechanical  .loads  that  the  system  must  withstand. 

Unsteady  operation  of  the  destructive  variety,  known  as  combustion 
instability,  was  first  encountered  in  1940.  At  that  time  a  British  group 
testing  a  small  solid-propellant  rocket  motor  observed  sudden  increases 
of  pressure  to  twice  the  expected  level,  enough  to  destroy  a  motor  of 
flight  weight.  Since  that  time  every  major  rocket  development  program 
has  been  plagued  by  combustion  instability  of  some  form.  These 
oscillations  in  the  combustion  chamber  can  have  several  detrimental  effects. 
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In  some  cases,  particularly  in  solid-propellant  rockets,  instability 
can  cause  the  steady-state  pressure  to  increase  to  a  point  at  which  the 
rocket  motor  will  explode.  In  liquid-propellant  rocket  chambers  experi¬ 
encing  unstable  combustion,  heat  transfer  rates  to  the  walls  considerably 
exceed  the  corresponding  steady  state  heat  transfer  rates,  resulting  in 
burn-out  of  the  walls.  If  the  chamber  can  survive  these  effects,  mechanical 
vibrations  in  the  rocket  system  can  cause  mechanical  failure  or  destroy  the 
effectiveness  of  the  delicate  control  and  guidance  systems. 

The  phenomenon  of  combustion  instability  depends  heavily  upon  the 
unsteady  behavior  of  the  combustion  process.  The  organized  oscillations  of 
the  gas  within  the  chamber  must  be  coupled  with  the  combustion  process  in 
such  a  way  as  to  form  a  feedback  loop.  In  this  manner  part  of  the  energy 
stored  in  the  propellants  becomes  available  to  drive  large  amplitude 
oscillations.  An  understanding  of  this  coupling  between  the  combustion 
process  and  the  wave  motion  is  necessary  in  order  to  predict  the  stability 
characteristics  of  rocket  engines. 

Combustion  instability  problems  in  liquid  propellant  rocket  motors 
usually  fall  into  one  of  three  categories  according  to  the  frequency  of 
oscillation.  Low  frequency  combustion  instability,  also  known  as  chugging, 
is  characterized  by  frequencies  ranging  from  ten  tc  several  hundred 
hertz,  nearly  spatially  uniform  properties,  and  coupling  with  the  feed 
system  of  the  rocket.  This  type  of  instability  is  less  detrimental  than 
other  forms,  and  the  means  of  preventing  it  are  well  understood.  Low  t 

frequency  instability  will  not  be  considered. 

A  second  type  of  combustion  instability,  which  is  less  frequently 
observed,  has  a  frequency  of  several  hundred  cycles  per  second.  This  j 
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type  of  oscillation  is  associated  with  the  appearance  of  entropy  waves 
inside  the  combustion  chamber. 

The  third  and  most  important  form  of  combustion  instability  is 
known  as  high  frequency  or  acoustic  instability.  As  the  name  suggests, 
this  type  of  instability  represents  the  case  of  forced  oscillations  of  the 
combustion  chamber  gases  which  are  driven  by  the  unsteady  combustion  process 
and  Interact  with  the  resonance  properties  of  the  combustor  geometry.  The 
observed  frequencies,  which  are  as  high  as  10,000  cycles  per  second,  are 
very  close  to  those  of  the  natural  acoustic  modes  of  a  closed-ended 
chamber  of  the  same  geometry  as  the  one  experiencing  unstable  combustion. 
High  frequency  combustion  instability  is  by  far  the  most  destructive  and 
is  the  type  to  be  considered  by  the  following  analysis. 

High  frequency  combustion  instability  can  resemble  any  of  the 
following  acoustic  modes:  (1)  longitudinal,  (2)  transverse,  and  (3) 
combined  longitudinal-transverse  modes.  Longitudinal  oscillations  are 
usually  observed  in  chambers  whose  length  to  diameter  ratio  is  much  greater 
than  one;  in  this  case  the  velocity  fluctuations  are  parallel  to  the  axis 
of  the  chamber  and  the  disturbances  depend  only  on  one  space  dimension. 

For  much  shorter  chambers  the  transverse  mode  of  instability  is  most 
frequently  observed.  Transverse  oscillations  in  rocket  motors  are 
^naracterized  by  a  component  of  the  velocity-perturbation  which  is 
perpendicular  to  the  axis  of  the  chamber  but  the  disturbances  can  depend 
upon  three  apace  dimensions.  Such  oscillations  can  take  either  of  two 
forms:  (1)  the  standing  form  in  which  the  nodal  surfaces  are  stationary 
and  (2)  the  spinning  form  in  which  the  nodal  surfaces  rotate  in  either  the 
clockwise  or  counterclockwise  direction.  Transverse  combustion  insta¬ 
bility,  particularly  that  resembling  the  first  tangential  mode,  has  been 
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frequently  encountered  in  modern  rocket  development  programs  and  has  been 
the  subject  of  much  current  research. 

Historic  Studies  In  the  Problems  of  Combustion  Instability 

Since  the  early  1950’s  much  experimental  and  analytical  research 
has  been  devoted  to  better  understanding  the  phenomenon  of  high  frequency 
combustion  instability.  Most  of  the  theories  presented  prior  to  1966  were 
restricted  to  circumstances  in  which  the  amplitudes  of  the  pressure 
oscillations  were  infinitesimally  small  in  the  linear  regime.  Prominent 
among  these  are  the  pioneering  studies  of  longitudinal  instability  by 
Crocco  [1]  as  well  as  the  studies  of  transverse  instability  by  Scala  [23, 
Reardon  [3],  and  Culick  [4].  A  complete  discussion  of  these  theories  is 
given  in  the  work  of  Zinn  [53  and  will  not  be  repeated  here. 

Although  linear  theories  provide  the  propulsion  engineer  with 
considerable  insight  into  the  problem,  their  applicability  and  usefulness 
in  design  is  limited.  The  linear  theories  cannot  provide  answers  to  such 
important  problems  as  the  limiting  value  of  the  pressure  amplitude 
attained  by  a  small  disturbance  in  the  case  of  a  linearly  unstable  engine, 
or  the  effect  of  a  finite-amplitude  disturbance  upon  the  behavior  of  a 
linearly  stable  engine.  In  the  latter  case  the  result  of  many  tests 
indicate  that  under  certain  conditions  the  introduction  of  sufficiently 
large  disturbances  into  a  linearly  stable  engine  can  trigger  combustion 
instability.  Another  shortcoming  of  linear  theories  is  the  fact  that 
their  predictions  cannot  be  compared  directly  with  available  experimental 
data;  for,  in  the  majority  of  cases,  the  experimental  data  is  obtained 
under  conditions  in  which  the  combustion  instability  is  fully  developed 
and  in  a  non-linear  regime.  Therefore,  theories  accounting  for  these 
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nonlinearities  associated  with  combustion  instability  are  needed.  A 
more  detailed  discussion  of  the  nonlinear  aspects  of  combustion  instability 
can  be  found  in  a  work  by  Zinn  [5]. 

In  the  field  of  finite  amplitude  (nonlinear)  combustion  instability, 
mathematical  difficulities  have  precluded  any  exact  solutions,  and 
approximate  methods  and  numerical  analysis  have  been  used  almost  exclusively 
for  this  reason  publications  in  this  field  are  scarce.  Notable  among  these 
is  the  work  of  Maslen  and  Moore  [6]  who  studied  the  behavior  of  finite 
amplitude  transverse  waves  in  a  circular  cylinder.  Their  major  conclusion 
was  that,  unlike  longitudinal  oscillations,  transverse  waves  do  not  steepen 
to  form  shock  waves.  Maslen  and  Moore,  however,  considered  only  fluid 
mechanical  effects;  they  did  not  consider  the  influences  of  the  combustion 
process,  the  steady  state  flow,  and  the  nozzle  which  are  so  important  in 
the  analysis  of  combustion  instability  problems.  Nevertheless,  pressure 
recordings  taken  from  engines  experiencing  transverse  instability  reveal 
the  presence  of  continuous  pressure  waves  similar  in  form  to  those 
predicted  by  Maslen  and  Moore. 

One  of  the  first  nonlinear  analyses  to  include  the  effects  of 
the  combustion  process  and  the  resulting  steady  state  flow  was  performed 
by  Priem  and  Guentert  [7].  In  this  investigation,  the  problem  was  made 
one-dimensional  by  considering  the  behavior  of  tangential  waves  traveling 
in  a  narrow  annular  combustor  of  a  liquid  propellant  rocket  motor.  They 
used  a  computer  to  solve  numerically  the  resulting  nonlinear  equations  for 
various  values  of  the  parameters  involved.  Due  to  the  many  assumptions 
involved  in  the  derivation  of  the  one-dimensional  equations,  the  results 
of  this  investigation  are  open  to  question. 
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The  successful  use  of  the  time-lag  concept  (see  Crocco  [1])  in  the 
lineal'  theories  prompted  a  number  of  researchers  to  apply  this  model  to 
the  analysis  of  non-linear  combustion  instability.  By  considering  a 
chamber  with  a  concentrated  combustion  zone  and  a  short  nozzle,  Sirignano 
[8]  demonstrated  the  existance  of  continuous,  finite-amplitude,  longitudinal 
periodic  waves.  These  solutions  were  shown  to  be  unstable,  however,  thus 
indicating  the  possibility  of  triggering  longitudinal  oscillations. 

Mitchell  [9l  extended  the  work  of  Sirignano  to  include  the  possibility  of 
discontinuous  solutions.  In  this  manner  he  was  able  to  show  that  the  final 
form  of  triggered  longitudinal  instability  consisted  of  shock  waves  moving 
back  and  forth  along  the  combustion  chamber.  Mitchell  also  considered  the 
more  realistic  case  of  distributed  combustion. 

In  the  analyses  of  Priem,  Sirignano,  and  Mitchell  the.  oscillations 
were  dependent  on  only  one  space  dimension.  One  of  the  first  researchers 
to  study  finite-amplitude  three-dimensional  combustion  oscillations  was 
Zinn  [5]  whose  work  is  an  extension  of  the  linear  transverse  theories  and 
the  analysis  of  Maslen  and  Moore.  Using  Crocco' s  time  lag  model  Zinn 
investigated  the  nonlinear  behavior  of  transverse  waves  in  a  chamber  with 
a  concentrated  combustion  zone  at  the  injector  end  and  an  arbitrary 
converging-diverging  nozzle  at  the  other  end.  In  this  case,  it  was 
necessary  to  extend  Crocco' s  burning  rate  expression  and  transverse  nozzle 
admittance  relation  to  obtain  the  appropriate  boundary  conditions  for  the 
case  when  the  flow  oscillations  are  of  finite  size.  As  a  result  of  this 
analysis  Zinn  was  able  to  prove  the  existence  of  three  dimensional 
finite-amplitude  continuous  waves  which  are  periodic  in  time.  In 
addition,  he  was  able  to  prove  the  possibility  of  triggering  combustion 
oscillations.  An  analytical  criterion  for  the  determination  of  the 
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stability  of  such  waves  was  derived,  but  because  of  its  complicated  form 
and  the  limited  capacity  of  available  computers  no  specific  numerical 
results  were  obtained. 

In  more  recent  years  other  investigators  such  as  Burstein  ClO] 
have  attempted  to  solve  numerically  the  equations  describing  instabilities 
that  depend  on  two  space  dimensions.  Although  the  resulting  solutions 
resemble  experimentally  observed  combustion  instability,  this  method 
requires  excessive  computer  time,  and  studies  of  this  type  for  three- 
dimensional  oscillations  will  have  to  await  the  development  of  a  much 
faster  breed  of  computers. 

In  a  recent  publication  by  Powell  [ll],  the  problem  of  analytically 
and  numerically  analyzing  multidimensional  non-linear  combustion  instability 
was  investigated.  The  problem  in  doing  this  is  that  a  system  of  non¬ 
linear  coupled  partial  differential  equations  whose  solutions  must 
satisfy  a  complicated  set  of  boundary  conditions  governs  the  phenomena  of 
combustion»instability.  These  boundary  conditions  may  describe  the 
unsteady  burning  process  of  the  wall  of  a  solid  propellant  rocket  motor; 
the  conditions  at  an  idealized  concentrated  combustion  zone  of  a  liquid- 
propellant  rocket  engine;  or  the  unsteady  flow  of  the  entrance  of  a 
converging-diverging  nozzle.  Previously,  in  an  effort  to  obtain  analytical 
solutions  to  various  combustion  instability  problems,  investigators  have 
been  forced  to  simplify  the  original  problem  to  such  an  extent  that  it  no 
longer  resembled  the  real  problem  that  originally  was  to  be  solved.  Powell 
proposed  a  method  to  perform  a  nonlinear  stability  analysis  with  relative 
ease.  This  methodf applicable  to  both  linear  and  non  linear  problems  with 
complicated  boundary  conditions,  was  a  modified  form  of  the  classical 
Galerkin  method.  The  Galerkin  method  [11]  is  an  approximate  mathematical 
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technique  which  has  been  successfully  employed  in  the  solution  of  various 
engineering  problems  in  the  field  of  acoustics.  Powell  used  this  method 
to  specifically  study  the  non-linear  behavior  of  combustion  driven 
oscillations  in  cylindrical  combustion  chambers  in  which  the  liquid 
propellants  are  injected  uniformly  across  the  injector  face  and  the 
combustion  process  is  distributed  throughout  the  combustion  chamber.  Based 
upon  the  results  of  his  second  and  third  order  theories,  the  following 
nonlinear  mechanisms  were  found  to  be  important  in  determining  the  non¬ 
linear  stability  characteristics  of  the  system:  (1)  the  transfer  of  energy 
between  modes,  (2)  the  self-coupling  of  a  mode  with  itself,  and  (3)  a  non¬ 
linear  combustion  mass  source.  Powell  found  that  the  self-coupling 
mechanism  was  important  in  the  initiation  of  triggered  instability,  while 
the  non-linear  driving  mechanism  was  important  in  the  determination  of  the 
final  amplitude  of  triggered  instability. 


Statement  of  the  Problem 

• 

In  this  thesis,  the  problem  of  velocity-sensitive  instability  will 
be  considered.  Based  upon  previous  work  on  this  problem,  only  transverse 
oscillations  will  be  considered  due  to  mathematical  simplicities.  Also, 
the  specific  geometry  of  the  combustion  chamber  to  be  analyzed  will  be 
annular  or  ring-like.  The  purpose  of  this  thesis  is  to  investigate  the 
mechanisms  which  cause  these  instabilities  due  to  the  combustion  process 


in  a  liquid  propellant  annular  combustion  chamber  and  attempt  to  state 
which  mechanisms  or  conditions  impose  the  greatest  effect  upon  stability 
of  combustion. 
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the  general  acoustic  wave  equation  for  non-linear  combustion  is  derived. 

In  this  derivation,  both  steady  state  and  deviations  from  the  steady-state 
conditions  are  considered  and  their  effects  incorporated  into  the  general 
acoustic  wave  equation. 

In  Chapter  3,  the  Galerkin  method  is  used  to  obtain,  from  the 
general  acoustic  equation  of  Chapter  2,  equations  governing  the  modal 
amplitudes  associated  with  the  first  two  modes  of  transverse  oscillation 
in  a  thin  annular  combustion  chamber.  These  equations  for  the  annular 
combustion  chamber  are  solved  numerically  by  the  use  of  a  Runge-Kutta 
program  for  various  conditions. 

In  Chapter  4,  a  set  of  approximate  equations  are  derived  from  the 
modal  amplitude  equations  presented  in  Chapter  3  by  use  of  the  two-variable 
perturbation  technique.  These  resulting  approximate  equations  are 
expressed  both  in  the  modal  amplitude  and  amplitude-phase  angle  form.  In 
this  chapter,  four  special  cases  are  presented  for  which  closed-form 
solutions  can  be  found.  These  four  cases  are  (1)  standing  wave — no 
combustion,  (2)  standing  wave — no  gas  dynamic  nonlinearities,  (3) 
traveling  wave — no  combustion,  and  (4)  traveling  wave — no  gas  dynamic 
nonlinearities.  For  problems  not  falling  within  the  above  categories, 
a  numerical  analysis  Is  employed  to  solve  approximate  equations. 

In  Chapter  5,  the  results  contained  in  the  previous  two  chapters 
are  discussed  and  compared.  Stability  limits  are  obtained  and  the  effect 
of  neglecting  various  physical  effects  are  discussed.  In  addition,  the 
accuracy  of  the  perturbation  method  is  evaluated.  A  summary  of  the 
research  contained  in  this  thesis  is  presented  in  this  chapter. 

In  Chapter  6,  a  statement  of  conclusions  is  made  along  with 
recommendations  for  future  research  in  this  area. 
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Chapter  2 


DERIVATION  OF  THE  GOVERNING  ACOUSTIC  WAVE  EQUATION 


In  order  to  investigate  the  non-l,inear  combustion  instabilities 
that  occur  in  liquid  propellant  rocket  engines,  one  must  start  with  the 
balance  laws  of  mass  and  momentum.  Also,  for  this  problem,  a  constitutive 
equation  was  formulated  relating  pressure  and  density.  Mathematically, 
these  principles  are  respectively 


5pft  ±  /*  •*.  _ . 
—3-  +  V  •  (p  u)  =  B* 


(2.1) 


/  *  A  ft  \ 


(2.2) 


ft  *2  ft 
P  3  a  p, 


(2,3) 


where 


p  -  gas  density 


t  -  time 


. >  g  j.  g  3  -> 

v  -  del  operator  of  the  system  -s^-  *  +  gy3F  3  +  gfTF  K 


u  -  velocity  of  the  gas 


B  -  fuel  drop  burning  rate  per  unit  volume 


p  -  pressure  of  the  gas 


a  -  constant  of  proportionality  (in  this  case  -  speed  of 
sound ) . 
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turner* 


The  *  representation  denotes  that  the  above  physical  quantities  are 
dimensional.  Equations  (2.1)  ~  (2.3)  are  based  on  the  assumption  that 
the  fuel  drops  serve  only  as  a  source  of  mass  for  the  gas  phase. 
Interphase  transfer  of  momentum  and  energy  are  neglected. 

Combining  equations  (2.2)  and  (2.3),  the  resulting  equation  is 

p(ft5T  +  u  ♦  $  u)  *  -  a2  $  p.  (2 


For  the  physical  situation  depicted  in  Figure  1 


Figure  1.  Schematic  of  a  Liquid  Propellant  Combustion  Chamber 
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A  convenient  non-dimens ionalizat ion  of  the  variables  is  as  follows: 


*  *  ,* 

P  *  PCP  (pq 


initial  density  of  gas) 


* 

-+ 

U  * 


*  + 
a  u 


ft 

?  s  T-? 

L 


*  *2  * 
p*  =  a  pQp 


*  * 

*  poa  n 

B  =  u  B* 


Substituting  these  non-dimensional  relations  into  equations  (2.1),  (2.3), 
and  (2.4),  the  results  are 


||  +  ?  •  <pu>  ■  B 


(2.5) 


P 


(2.6) 


pap 


(2.7) 


where  the  unstarred  quantities  are  dimensionless. 
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Dividing  through  by  density  p,  equation  (2.6)  becomes 


3u  .  •+  *  vp 

Tt  +  u  *  *  u  s  • 


(2.8) 


Since , 


~  =  ^  6t  P, 


the  governing  equations  can  be  summarized  as 


||  t  ?  •  (ou)  =  B 


(2.9) 


«s.  -ttap 


(2.10) 


p  =  P. 


(2.11) 


It  will  now  be  shown  that  to  the  order  of  approximation  inherent 


in  these  equations,  the  flow  is  irrotational ,  that  is  5  x  u  =  0.  To  do 


this,  take  the  curl  of  equation  (2.10)  and  set  it  equal  to  zero.  The 


resulting  equation  becomes 


»*($♦«•  *«)■-***&.  *0. 


(2.12) 


Since  the  curl  of  any  gradient  is  zero.  This  may  be  rewritten  as 


?  x  *  (i  •  ?  i>  =  o. 


i.  1  - -  - 


(2.13) 


Therefore , 

^  x  (u  •  ^  u)  =  7  x  [^(*5u2)  ■  u  x  2]  ,  ,  (2.17) 

Recognizing  that  the  curl  of  any  gradient  is  zero,  equation  (2.17) 
reduces  to 

$  x  (u  •  ^  u)  =  -  7  x  (u  x  j^)  .  (2.18) 

Using  the  vector  identity 


v  ,  (1 «  I)  ■  (I  .  W  -  S<?  -  h  -  a  ■  hi  *  W  ■  I) 
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equation  (2.18)  can  be  expressed  as 

»  *  <»  .  *  «>  «  -  t<s  •  *>  s  -  s<*  •  u>  -  <s  •  ?>a 

t  u  (5  -3)3  .  (2.19) 

Therefore,  equation  (2.13)  becomes 


!£  -(M)St  3(V  .  u)  +  (u  •  8)8  -  8(8  •  8)  3  0.  (2.20) 

Equation  (2.20)  can  now  be  modified  by  using  the  definition  for  the  total 
(comoving)  derivative  which  is 


D£  dfi 
Dt  '  3t 


+  u 


8). 


Substituting  this  expression  into  equation  (2.20)  and  simplifying,  the 
resulting  equation  becomes 


j|.  8  •(?!)  -  (8  $>.«♦«<»•  8). 


(2.21) 


Rewriting  u  ($  •  8)  as  u  [$  •  ($  x  u)]  which  is  zero  since  th«  divergence 
of  the  curl  of  any  vector  is  zero,  equation  (2.21)  becomes 


§f .  a .  c?  -  (8 1>  ■  : . 


(2.22) 


The  implications  of  this  equation  for  a  fluid  starting  from  rest  are  as 
follows.  At  the  initial  instant  of  time  (t  3  0),  the  vorticity  of  any 


fluid  particle  will  be  zero.  Thus,  the  time  derivative  of  the  vorticity 
of  the  particle  will  be  zero,  implying  that  »  0  at  t  *  0.  Since 

•to 

H  ■  0  and  jyj»  *  0  at  t  *  0,  it  follows  that  ?!  *  0  at  the  next  instant  of 
time.  By  induction,  it  can  be  shown  that  5$  =  0  for  all  time  unless  the 
velocity  gradient  becomes  infinite  for  any  t  =  0.  It  is  assumed  in  what 
follows  that  this  does  not  occur  and  the  flow  is  treated  as  irrotational. 
Since  irrotationality  has  been  proven,  the  velocity  vector  u  can 
be  expressed  as 

U  =  $  ip  (2.23, 

where  1 1>  is  the  velocity  potential.  Substituting  equation  (2.23)  into  the 
left  hand  side  of  equation  (2.10),  the  result  is 

f§  ♦  u  •  ?  I  .  |f  ♦  ?Ow2)  -  I  X  5 

=  (5  *)  +  *[<,<?  *)^|  -  ?  *  X  8  .  (2.24) 

For  Irrotational  flow  ( *  0),  the  right  hand  side  of  equation  (2.24) 
becomes 


•  ?*j)  _ 

Therefore,  equation  (2.10)  can  be  written  as 


(2.25) 


^ftt-  +  *s(Vij>  •  $(j/)  +  In  pj  =  o  . 


(2.26) 


1? 


Spatially  intagrating  aquation  (2,26)  produces 

i 

+  *9  h  '  'h  ♦  *>t  p  ■  o(t)  (2,27) 

where  «(t)  is  a  function  of  integration.  From  equation  (2.23),  it  can 
be  seen  that  an  arbitrary  function  of  time  can  be  added  to  without 
affecting  the  result  for  u.  Thus,  a(t)  could  be  absorbed  into  The 
sam«  thing  is  accomplished  by  setting  a  *  0  which  results  in 


*  ’  -  If  -  (2.28 

or 

p  »  e  '  '  .  (2.29 

Thus,  p  and  u  are  both  known  as  functions  of  4*.  From  equation  (2.9),  the 
governing  equation  for  4<  can  be  written  symbolically  as 


||  +  t  •  %  *  B  (2. 30. a) 

-fe  ♦  *  h  *  *♦) 

p  ■  p  *  e  'l  ‘  ,  (2. 30. b) 


Rather  than  combining  these  quantities  immediately,  it  is  convenient  to 
first  make  further  simplifications  based  on  the  nature  of  the  physical 
problem  that  it  is  desired  to  analyse. 


18 


Steady  State  Solution 

First,  the  steady  state  solution  of  equations  (2,30)  corresponding 
to  purely  axial  motion  will  be  found.  Define  the  steady-state  velocity 
potential  $  by 

$  *  c$(a)  (2.31) 

where  c  (assumed  small)  is  the  measure  of  the  deviation  of  the  density 
from  its  initial  value  (see  equation  2.32  below).  The  bar  notation  will 
represent  steady-state  conditions.  The  steady-state  burning  rate  w  is 
defined  from 


B  =  co(e).  (2.32) 

While  many  other  situations  are  possible,  attention  will  be  confined  in 
the  present  work  to  the  case  when  to  =  0(e).  To  indicate  this  let 

C  =  e3  (3  =  0(1)).  (2.33) 

ThU3,  the  burning  rate  B  can  be  expressed  as 

B  =  to.  (2.34) 

Equation  (2.30.b)  can  now  be  written 


P 


(2.35) 


Using  the  Taylor  series  expansion  for  the  exponential  function  and 
retaining  only  the  first  two  terms,  equation  (2.35)  becomes 


f  *  1  -  •}*' 


(2.36) 


Substituting  aquations  (2.31),  (2.34),  and  (2.36)  into  aquation  (2. 30. a) 
and  dividing  the  result  by  c  yields 


(2.37) 


£jj.  _  2  fefjY  - 

a??  f c  [z tj  \&£r  *  •  * 


(2.38) 


Retaining  only  terms  of  0(1)  produces 


£i  *  3. 

dz<- 


(2.39) 


For  simplicity,  only  the  case  of  uniformly  distributed  combustion  (i.e. 
o  -  constant)  will  be  considered.  Thu3,  integrating  equation  (2.39)  one 
obtains 


“2  so  2  +  C 
dz  1 


(2.40) 


where  =  u  is  the  steady  state  velocity  of  the  gas. 


At  the  injector  (Z  *  0),  u  *  0.  Thus,  Cj  *  0  and 


u  *  4^  =  OS. 
dz 


(2.41) 


pimp . 


\ff  ytrorpl*. stj.  i "  )1’TT 
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Deviations  from  Steady  State 

It  la  now  desired  to  investigate  the  stability  of  the  steady 
state  solution  discussed  above.  Toward  this  end*  an  additional 
velocity  potential  related  to  perturbations  from  the  steady  state  is 
defined  by  the  equation 

♦  *  e  C5  +  $(*,  y,  *»  t)3. 

A  perturbation  burning  rate  B  is  also  defined  by  the  equation 

B  s  w  +  ew. 


(2.42) 


(2.43) 


It  is  assumed  that  w  *  0(t)  and  this  is  indicated  by  defining  a  function 
o  such  that  o  =  0(1)  and  w  *  oe.  Then  equation  (2.43)  becomes 


B  =  e(o  +  eo)  , 


Taking  the  gradient  of  equation  (2.42),  one  obtains 


(2.44) 


Vi|>  =  e[$  $  +  $  $3 


(2.45) 


or 


i>  *  e[u  a„  +  $  $3 . 

m 


From  equation  (2.42),  the  time  derivative  of  iji  can  be  expressed  as 


(2.46) 


3£  B  e  3 

at 


3* 

at  . 


(2.47) 


Substituting  th«  aquations  (2.46)  and  (2.47)  into  aquation  (2.30.b) 
and  simplifying,  ona  obtains 
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p  s  p  *  q 


(2.48) 


3 

Expanding  (2.48)  in  a  Taylor  series  and  neglecting  terms  of  0(e  )  and 
higher  produces  the  expression 


p  =  p  =  i_eli+e2  ^(G2  ■* 


(2.49) 


Substituting  equations  (2.42),  (2.44),  and  (2.48)  into  equation  (2. 30. a) 
and  dividing  the  result  by  e  leads  to 


+  £ 


+  . 


+  .  $*)  »g  |4-+  *s  I-  (I1!2! 

9  97  3z3t  ’  3t  \3tj  J 

6  Tt  +  e2  +  ^  ^  "  u 


(-  •*  it 

+  e2V(-is(u2  +  $$  •  $$)j-  u  +  h  J  *  5  +  eo  . 


(2.50) 


Neglecting  all  terms  of  0(e2)  and  higher  and  recalling  from  the  steady- 
state  solution  that  u  *  =  oz  and  ^  *  o  yields 
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0  -  v2*  +  c  [ja  <$♦  •  ^*)  +  u  }if-  -  |1  |i|.  +  |i  ^2^ 

+  It  5  *  G  lilt  +  (**  *  *  It)]  *  -oc*  (2-51) 

Substituting 

*  ~  ($♦  •  ft)  *  *♦  *  $  H*  <2.52> 

into  equation  (2. SI),  results  in 


|^  -  72^  +  e  [2^  *  ^  ff)  ♦  2  «  fift  +  ff  5 


a  -oe 


(2.53) 


where  only  terms  of  0(1)  and  0(c)  have  been  retained.  Equation  (2.53)  can 

j2i 

be  further  simplified  by  observing  that  72<f>  =  +  0(c). 

Thus,  the  last  term  of  equation  (2.53)  can  be  written 


Siflii 

at  \3tT 


+  0(e)  - 


0(c2 ) . 


Since  the  other  terms  of  0(e2)  have  already  been  neglected,  consistency 
requires  that  this  term  be  deleted  and  the  equation  be  rewritten  as 

ft*  “  72*  +  6  ^  (*♦  *  ^  It)  +  2G  lift  +  ft 


5  = 


-oc 


(2.54) 
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In  this  thesis*  attention  will  be  confined  to  transverse  instability. 

For  this  situation 


♦  ■  ♦(*,  y*  t).  (2.55) 

Therefore,  equation  (2.54)  becomes 


(2.56) 


To  account  approximately  for  frequency  changes  due  to  baffles,  nozzle 
shapes,  etc.,  a  correction  term  of  the  form 


6  K*2(!t£)  (2.57) 

was  introduced  into  equation  (2.56).  This  form,  one  of  many  possible,  was 
chosen  so  that  the  linearized  form  of  equation  (2.56)  would  reduce  to  Love's 
equation  for  a  one-dimensional  problem.  This  linearized  form  of  (2.56)  is 

3t7  "  3x*  '  6  Kax2at 7  -  °*  (2.58) 


Thus,  it  can  be  seen  that  the  value  K will  affect  the  acoustic  frequencies. 
Physically,  this  is  the  purpose  of  baffles,  nozzle  shapes,  and  other 
physical  parts  of  the  combustion  chamber.  Therefore,  inserting  the 
correction  term  into  equation  (2,56),  the  resulting  equation  becomes 
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where  K  is  the  correction  factor.  This  non-linear  wave  o qua t ion  will  be 
the  basis  for  numerically  and  analytically  investigating  the  transverse 
combustion  stability  problems  occurring  in  liquid  propellant  rocket 
engines . 


) 
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Chapter  3 

DERIVATION  OF  WAVE  EQUATIONS  BASED  UPON  AN 
ANNULAR  COMBUSTION  CHAMBER 

In  Chapter  2,  there  were  no  restrictions  concerning  the  geometry 
of  the  combustion  chamber  in  the  derivation  of  the  acoustic  wave 
equation.  In  this  chapter,  however,  a  set  of  equations  will  be  developed 
based  upon  a  narrow  annular  combustion  chamber.  A  typical  cross-section 
for  such  a  combustion  chamber  is  shown  in  Figure  2  below  in  dimensional 
and  dimensionless  form. 


(a)  Dimensional  (b)  Dimensionless 

Figure  2,  Dimensional  and  Dimensionless  Form  of  a  Circular 
Cylindrical  Combustion  Chamber 


In  Figure  2  (a),  the  dimensional  quantities  are 

ft 

r  -  radius  of  a  typical  point  in  the  combustion  chamber 
ft 

R  -  inside  radius  of  the  combustion  chamber 
ft 

b  -  thickness  of  combustion  chamber's  cross-section. 
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In  Figure  2  (b),  the  dimensionless  quantities  are 

r  -  non-dimensional  radius  of  a  typical  point 
A 

R  -  |=  1 
R 


The  first  major  assumption  to  be  made  in  the  geometry  of  the  combustion 
chamber  is 


b 

R 


<< 


1 


(3.1) 


which  states  that  the  circular  cylinder  can  be  thought  of  as  a  thin 
(ring-like)  annulus. 

Define  the  characteristic  length  L*  by 

L*  =  R*.  (3.2) 


In  restricting  the  analysis  to  an  annulus,  a  transformation  to  polar 
coordinates  is  convenient.  Recall  that  the  gradient  and  Laplacian 
operators  in  polar  coordinates  are 


iL. 

3r 


r  30 


1_ 

3z 


.32  13  .  1  32  .  32 

3r^  r  3r  r^"  302'  3z^ 


(3.3) 


The  second  major  assumption  for  the  simplification  of  the  velocity 
potential  is  restricting 
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♦  a  *(e,  t) 


r  *  1  . 


(3.4) 


Therefor®,  using  the  operators  of  equations  (3.3)  on  the  function  of 
equation  (3.4),  the  results  are 


h 


■*  34 

e8  30 


(3.5) 


Substituting  the  results  of  equation  (3.5)  into  the  general  acoustic  wave 
equation  (2.58),  the  modified  wave  equation  becomes 


+  e 


o  3^4 

30  '  303t 


K34*  1  - 

*3^36^1  * 


oe  . 


(3.6) 


Now,  express  the  velocity  vector 


u  =  u  +  u 


(3.7) 


where  u  -  steady-state  velocity  vector 
u  -  perturbation  velocity  vector. 

From  the  steady  state  solution  in  Chapter  2,  the  velocity  vector  was 
defined  as 


Q  *  c  [§!  *z]  . 


(3.8) 
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Define  the  perturbation  velocity  vector  by 


-n 

u 


(3.9) 


Substituting  equation  (3.8)  and  (3.9)  into  equation  (3.7)  and  using 
equation  (2.23)  results  in 


u  =  e 


i 

dz  “z 


+  e 


3$ 

Te 


ee  -  7*. 


(3.10) 


To  determine  only  the  transverse  velocity  component  of  the  perturbation 
velocity  vector,  subtract  the  perturbed  velocity  component  along  the 
axial  (z)  direction  of  the  chamber  from  the  total  perturbation  velocity 
vector.  Thus, 


-  u  e  . 
z  z 


(3.11) 


In  this  case,  since  u  *  u(0,  t)  only,  there  is  no  perturbed  velocity 
component  in  the  axial  direction;  therefore. 


Ut'  =  e  38  eS. 


3<P 


(3.12) 


It  is  now  desired  to  find  the  burning  rate  a  in  terms  of  the  parameters 
in  the  wave  equation.  To  obtain  this  expression,  assume  velocity  sensitive 
combustion  with  no  history  effects.  Mathematically,  the  burning-rate 
function  for  velocity-sensitive  combustion  will  be  expressed  by  the  purely 
phenomenological  equation 
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o  s  u  n  f 


(3.13) 


where  n  is  called  the  interaction  index. 

Using  the  derived  results  for  the  general  time-delay  integral 
(discussed  in  Appendix  A),  the  burning  rate  with  history  effects 
accounted  for  by  a  simple  time  delay  is 


o  *  u  n 


(3.14) 


where  the  subscript  t  represents  the  time  delay.  For  simplicity,  it 
will  be  assumed  that 


Then,  the  burning  rate  can  be  expressed  as 


(3.16) 


where  j  =  0  -  no  time  delay 

1  -  time  delay  . 

Therefore,  substituting  equation  (3.16)  into  equation  (3.6),  equation 
(3.6)  can  be  rewritten 
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+  nw 


(3.17) 


There  is  no  closed  form  solution  of  equation  (3.17)  that  appears  likely. 
The  main  purpose  of  the  present  work  is  to  determine  the  modifications  of 
solutions  of  the  usual  acoustic  wave  equations  that  are  caused  by  the 
presense  of  the  nonlinear  terms  multiplied  by  e  in  equation  (3.17). 

Thus,  rather  than  attempt  a  finite  difference  numerical  solution  of 
equation  (3.17),  the  following  procedure  was  adopted. 

The  solution  is  represented  by  the  Courier  series 

$(9,  t)  =  f ^ ( t )  cos  9  +  f2(t)  cos  26  +  g^(t)  sin  6 

+  g2(t)  sin  20  +  .  .  .  (3.18) 


and  initial  conditions  are  chosen  such  that  in  the  absence  of  the  nonlinear 
terms,  the  exact  solution  can  be  formed  using  only  the  first  two  terms  of 
the  Fourier  series.  Because  of  the  quadratic  nature  of  the  non-linearities, 
the  second  two  terms  in  equation  (3.18)  represent  a  complete  first  order 
correction  to  the  acoustic  solution  due  to  non-linear  gas-dynamic  and 
combustion  effects.  Only  the  first  four  terms  in  equation  (3.18)  are, 
therefore,  retained  and  the  approximate  solution  determined  by  this  method 
is  the  simplest  one  capable  of  illustrating  the  influence  of  the  nonlinear 


terms.  The  approximation  can,  of  course,  ba  improvad  by  rataining 
additional  tarma  in  aquation  (3.18)  but  this  ia  not  invaatigatad. 

Substituting  aquation  (3.13)  into  aquation  (3.17)  and  using 


tha  multiple  angle  formulas  to  simplify  tarma  containing  products  of 
trignomatric  functions,  ona  obtains 

{d2f  _  df  r  df  df  dg  dg  1 

4  fi  +  “  St1  *  2c  [f2  at1  *  fi  Hi*  *  *2  371  *  *i  3t  J 

*  *  dt^  +  2nc“  [flf2  *  *182]  -  2^"“  [fl,f2t  f  glT82]|C08  9 


* 

d2g  dg  r  df  dg  df  dg  T 

■nitT1  +  «1  +  “  dt1  +  2e[g2  dt1  *  fl  dt*  "  gl  dt*  “  f2  dtl 

+  KC  ^  t  2  neSjf^  -  f2gA]  -  2jein  [f^  -  f2TglT]|  sin  0 


♦for*  +  4  f 2  + « dt^ +  e 


'  dg,  df,l 

?i  dT1 "  fi  dtM  *  4Kc 


t  VJen  [gL2  -  fj2]  -  Wien  -  fj2] 


cos  2G 


+fe?*  +  4g2  +  “  "  C  gl 


-r~i"  +  f 


dg  d2g 

1  dt1  +  4Ke  dt2 

m 


*  iinE  [fi 


sin  20  +  .  .  ,»0. 


(3.19) 


Equation  (3.19)  ia  a  summation  of  terms  composed  of  some  function  of 
time  t  and  a  term  containing  0  variation.  Since  the  equation  must  be 
valid  for  all  values  of  0,  each  of  the  time  dependent  coefficients 
of  the  0-terms  must  individually  be  equal  to  aero.  Therefore,  four 
ordinary  differential  equations  governing  the  time-dependent  modal 
amplitudes  and  emmergo  from  this  analysis  as  the  governing 

equations  to  be  used  for  analysis  of  instability  in  an  annular  combustion 
chamber.  These  equations  are 


!f  -  df,  f  dfi  df2  dgi  dg21 

-?1  +  fi  +  “  dt1  +  2e  [f2  dt"”  +  f  1  dt~  +  g2  dt~  +  81  dFj 


**  dt^  +  2nEtu  [flf2  +  filg2]  “  2;)enui  [flTf2T  +  ^*1  "V 82t}  =  0 


(3. 20. a) 


d2g 


dg  r  df 

+  gi  +  “  3t  *  2cLg2  dt1  +  h 


dg  df  dg 

dt  gl  dt  r2  dt 


.  Kc  .  2nd  [fltj  -  f,*,]  -  2jcin  [flTgJt  -  fjtsJ  *  0 

(3.20.1 


St*2,  +  4f2  +  “  dt^  +  6  h  dt**  ‘  f!  dt1  +  4KE  dt? 


t  Hein  J^2  -  f12j  -  ^jeun  Jgj2  -  f^j  *  o 


(3.20.c) 
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d2g  ^  dg  df  dg  d2g 

dt*  *  4*2  +  5  dt2"  *  c  *1  dt1  +  fl  dt1  +  4Ke  df?" 


[fi8i]  ♦  J'5*  ['uh'\  '  0  • 


(3 .20 .d) 


In  the  following  work  only  instantaneous  combustion  will  be  considered. 
Thus,  the  appropriate  equations  are  equations  (3.20)  with  j  *  0.  These 
equations  are  recapitulated  below. 

d2f  df  df  df  dg  dg 

3^  *  fi  *  “  371-  *  2c  f2  at1  *  fj  at*  *  h  dt1  +  gi  at* 

_  •> 


d2f  r  -| 

Kc  +  2ncw  +  gjgJ  *  0 


(3. 21. a) 


d2g  _  dg  df  dg  df  dg 

dt^ +  gi +  “  2?1  +  2e  g2  at1  +  fi  at2  ■  gi  at2  ’  f2  at1* 


+  Ke  +  2neoj  If^  -  f2g 


.]■ 


(3.21.b) 


d2f 

+  4f2 


r 

df 

+  “  dt2,  +  c  [fi 


dg  df 

l  _  f  _ t 

dt  rl  dt 


d2f  r  T 

4Ke  dt*2,  +  *sc“n  gi2 "  fij  s  0 


(3.21.c) 


aap.VTOiimiyv'flVil  “ 


',vwrw.  '.vtr  ■»TT^,wvi  ■ 


*7,pw<tT’MTn •'•«’ "s^rr n™:  !'"'•“  ‘ 


f*  m 

d2g  dg  df  dg 

jti2,  +  "«j  +  “  at  •  *  *1  at  *  fi  at^ 

»*  *■ 


d2a,  _  r  -j 

4Ke  ^y2-  -  umc  UjfJ  *  0 


(3.2i,d) 


The  equations  of  (3.21)  were  solved  numerically  by  the  use  of  the 
quartic  (fourth-order)  Runge-Kutta  method.  To  use  this  method,  the 
equations  of  (3.21)  are  modified  by  defining  the  quantities 


dt  *  al 


^tLs 


dT*hi 


dg, 

— — 2*  S  h 

dt  2  . 


(3.22) 


Substituting  these  expressions  into  equations  (3.20)  and  solving  these 
equations  for  the  highest  derivative  (in  this  case  -  second  order),  we  get 


dt1  *  {ffl  “  “^al^  ”  2e±(f-'2Cai >  +  f1(a2)  +  g^bj) 
+  g1<b2))  -  2neu(f1f2  +  gjg^J  /(I  +  Ke) 


1  ■  uwfiPnspwwwfH 


<us  r  / 

dt1  *  f*l  *  "(V  •  2cl\g2(al)  +  W  -  8l(a2} 

-  fjCb*))  "  Znewtfjgj  -  f2gl^]  +  **  * 


da 

_ 2 

dt 


£-4f2  -  u(a2)  -  «i(g1(b1)  -  f^a^J 


-Vwnjgj2  - 


/(I  +  4  Ke  ) 


db 

m 

dt 


*  [~4g2  -  w(b2)  +  ei|g1(a1)  +  fjCbjjJ 
+  une  (f181)j  /(I  +  4  fc) 


35 


(3.23) 


where  i  is  the  gas-dynamic  index. 

By  the  development  of  a  computer  program  incorporating  the  Runge- 
Kutta  algorithm  which  can  solve  systems  of  first-order  ordinary  differen¬ 
tial  equations,  the  eight  equations  (3.22)  and  (3.23)  were  numerically 
solved  for  the  eight  variables  a^,  a2,  bj,  b2,  f^  f2,  and  g2. 
Different  cases  involving  varying  the  gas-dynamic  index,  interaction 
index,  the  correction  variable  (K),  and  the  order  term  (epsilon)  will  be 
discussed  and  compared  with  the  perturbation  method  of  solution  in  a 
later  chapter.  In  Appendix  B,  a  sample  program  listing  this  calculation 
appears . 


Chapter  4 

TWO-VARIABLE  PERTURBATION  METHOD  APPLIED  TO  THE 
ACOUSTIC  WAVE  EQUATIONS 

In  this  chapter,  a  set  of  approximate  equations  will  be  developed  ( 
from  the  governing  equations  for,  the  modal  amplitudes  (3,21),  by  the  use 
of  the  two-variable  perturbation  method.  The  two-variable  method  is  well 
suited  to  this  type  problem  since  one  expects  the  solution  to  consist  of 
sinusoidal  functions  with  slowly  varying  amplitude.  Applying  this  method, 
define  two  variables  representing  time 

5  »  t 

n  a  6t  .  (4.1) 

Therefore,  the  four  modal  amplitudes  would  now  be 
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O 


G 


k. 


! 

j 


O 


o 


j-jS,  a  t  2( 

dt2  B52 


(4.4) 


where  Z  *  f^,  fj,  g^,  g2  respectively  for  each  of  the  above  equation#. 
By  substituting  equations  (4.3)  and  (4.4)  for  eaoh  modal  amplitude  into 
equations  (3.21)  and  keeping  terms  only  of  0(1)  and  0(c),  the  resulting 
equations  become 

a2*,  a2f ,  af,  af,  af,  ag, 

TT  *  'x  *  «P  *•  TC  *  "»  ar  +  8,1  W  *  ***  aT 

ag2  a^f, 

+  2gi  a T  +  K  aF"  +  2nw(fif2  +  0 


a2g,  a2g,  ag,  af,  ag5  af2 

ir*1!*  e[2H^  ”>r  *  ir  *  !fi  ir '  2gi  sr 


3gi  3iigi 

2f2  alp  +  K  “35”  +  2n<*»(fi62 


f2gl)]  “  0 


a2f~  a2f~  af2  ag, 

W  +  4f2  +  et2  Ti an  +  0  ‘aT*  +  8i  aT* 


a2f. 


♦  4K  jg-  +  h  un(g12-f12)]  =  0 


a2g2  a*g2  ag2  af. 

at2”  *  4g2  +  eC2a?^  +  °  —  ’ 


-  f. 


ag. 


as  ®i  as  1  ae 


+  4K- 


a2g,  _ 

:^2~  -  u»n(f1g1}]  =  0  . 


(4.5) 
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From  the  straight-forward  perturbation  method,  define  the  modal  amplitudes 
by  the  series  expansions 

fl  “  fl0(c*n)  +  c  fll(c»n)  +  •  •  • 

f2  *  f20^*n>  +  c  f2i^»n)  +  .  .  . 

8i  8  gio*5»n*  +  c  8nU»n)  +  •  •  . 

~  g2o^»n)  t  6  t  «  »  »  »  (4.6 ) 

Again  by  applying  the  rules  of  differentiation,  it  can  be  shown  that 

3Z  3T  ,  3K 
35  35  6  35 

liZ  _  32T  32K 

ae2  '  35"  6  W 


32Z  .  3^T_  32K 

353n  “  353n  c  353n 


(4.7) 


where 


Z  —  f^»  ^2*  gl*  ®2 


and 


^  —  f  i  f 2^»  8n *  g21^  respectively* 


Substituting  the  expressions  of  (4.6)  and  (4.7)  into  equations  (4.5)  and 
keeping  terms  only  of  0(1)  and  0(e),  the  resulting  equations  become 


■»^"nT?*rv jti 


** 'io  .  -  "10 


W~  *  fio  +  cC5?r~  *  fu  +  2  "S^TT  *  ®  IT"  + 


'2o  as 


+  2f  Ilsi  +  2a  lfl2.  ♦  2ff  !fao  +  k  a2fl° 

2fio  as  +  2g20  as  *  2gio  as  *  K 


+  2nw(f1Qf20  +  g10g20)3  *  0 


*2g10  a2gll  a2g10  -  3*io  Wla 

as7"” +  gio  +  cCIf*“~  *  gn  ♦  2  TO”" *  0  IT"  *  2g2o 


+  2f  3g20  -  5»  3f20  2f  3g10  K  3Zg10 

+  2fio  “  2gio  TT~  “  2f20  ~  +  K  W~~ 


+  2nw(f10g20  -  f20g1Q)]  *  o 


l!*2U4f  ,erl!fj!>t4f  ,  -  aif20  -»20 

IT2  4f2o  *  c  as*  +  4f2i  2  H"a"n  +  as 


sio  as 


»a.  .  ax  - -  ^ 

■io  Is  +  4K  W22’  *  fn(gio2  “  fio2)] 


32g?0  32g21  32g20  -  3g2n 

asz  *  4g20  *  ^Is2  *  ,+s2i  +  2  a s an  +  0  “aT* 


3^10  3g10  32®20  — 

‘gio  as  fio  as  +  4K  ~&sz  “  wn(fiogio)-!  s  0 


(4.8) 


By  separating  the  terms  of  0(1)  and  0(e)  in  the  equations  of  (4,8)  and 
equating  both  seta  of  tern)3  equal  to  zero,  the  resulting  equations  become 


a2fio 

a?2  +  fio  *  0 


ic*-  ♦  «!o  -  o 


+4f  *  o 

20 


3**20 

HS—  *  0 


f  =  -jlffiS.  _^io  sf10  Jfj0 
T?~  u  *5iSr  »Mlr*  -  2fl0_2£ 


(4. 9. a) 


2_  3gI0 
*2820Tr  28 


3g20  32fio 

10 1?-  -  *  ijr-  -  *■»*„»» ♦  g10g20) 


iiu  ♦  g  .  -  2  i!fi«  -  9gio 

S11  '  -  2  -  0  tf-  -  Jg; 


1-10  _  MJ0 

20  3£  2fio  a  f“ 


"81°  ^  +  2,io  TT1  -  K  5^  -  2"Sf10«20  -  f20gio) 


9  f2l  , ,  32*V-  3*f„ft  aa 

sp-t4f2i'-2i75r-»— -Sa 


*10  K  10  i7 


-"K  -  ^h(g102  .  fio2, 


^♦4g., ,  .jiffs  -*20 

21  *  ac3n  a  +  * 


10  35 


.  -  «10  3  Soq 

10  ~  '  4K  W~  *  n“uio*io) . 


(4.9.b) 


ipW.IJl|lllW|^W«!WIWW>pillWlli*l|l',lipiWpiilliWIHIIlll»l''>'-.-11'1  lu'1'  q»^SiT-H»^»»wrTnwy«^  .™,i.,„r,.-^,,-^,r„™t  -,,,,, 
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The  equations  of  (4. 9. a)  are  linear  second-order  differential  equations. 
Therefore,  it  can  be  shown  that  assuming  the  appropriate  form  of  a  solu¬ 
tion,  the  results  become 


f 

10 

If 

> 

H* 

w 

cos 

K  + 

B^n) 

sin  5 

*10 

1  Vn) 

cos 

5  + 

B2(n) 

sin  5 

f20 

*  A3(n) 

cos 

25  + 

B3(n) 

sin  25 

g20 

=  A4(n) 

cos 

25  •*• 

B4(n) 

sin  25  , 

(4.10) 

Substituting  (4.10)  into  (4.9.b)  and  using  the  multiple-angle  formulas 
yields 

^2^11  f  ?^1  1  —  1 

“  +  fii"“T  ST  -  7  V  +  KA3  ♦  bib3> 

-(B^Bg  +  A^)  +  -|(A2A4  +  B2B4)  -  (A2A4  +  B2E>4) 

-  |KB1  +  nuC^A^Bg  -  AgB1)  +  ^A^  -  B^)]]  sin  5 
!"dBl  1-  1 

“2[dTr  +  2«B1  +  *A3B1  -  AlV  +  <A1B3  -  A3Br) 

+  i  (A4B2  -  A2B4)  +  (A2B4  -  B2A4)  -  iKAx 

tnojC^A^Ag  +  B^Bj)  +  y  ^A2A4  +  B2B4^1  cos  £  ♦  •  •  • 


» 

♦  gu  ■  -2  - 

Ml 


S7-|”‘!’  ?<A1A4  *  W 


“(A1A4  +  B1B4^  +  ^A2A3  +  ®2B3^  ”  "2^A2^3  +  B2B3^ 


-  i  KB2  +  nwC^A^B^  -  A^)  -  ^(AjBg  -  AgB^jj  sin  5 


K  x— 

“2[d^  +  ¥*2  * 


¥\*1  ~  Vl>  +  (A1B4  '  W 


-(A2B3-B2A3)-|<A3B2-A2B3)-fKA2 


nuTt^A^  +  BxB4)  -  -|<A3A2  +  B2B3)]j  cos  5  +  .  .  . 


+  4f21  =  -2^-: 


dA3  1- 
ldrT  ’  Tj(2A3) 


+  ~  t|<B22  -  A22)  -  -  Ax2)3  +  2K(-4B3) 

♦  in»[A2B2  -  AJS^sin  25  -2^  +  ^(ajJt^AjBj-A^) 


COS  25  +  .  . 


+2K(-4A3>  +  ~nw[-|<A22  -  B22)  -  i<A^2-B1 2)]j  cos  25  +  .  . 

^T1  *  “«21  -  -{-4s1  *  -  &<B182  -  A1A2^ 

* 

+  ¥*1*2  *•  AlA2)]  t  2KC-4B4]  -  ^n^A^+A^J  sin  25 

* 

"  dB 

-2  J3n1  *  ?*2V  - 


2^2^B2A1+A2B1^+2^B2A1+A2B1>  ^ 


i 


WT' 


WPHPPjWPBiW  J'rj^Tp 


r  ,TR»W  .MHVTWT'V;’- 


rJ^.TT»T7fwji,»»i»-Tf^B,"*v,w  wirtjrn ' ;;  Tyr.ir*-  :,‘T'r! 


C- 


c 


♦  2KC-4Alrl  -  ■jnuiC^A 


1^1*2  -  WJ] 


COS  2£  +  .  .  . 
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C4.ll) 


where  +  .  .  .  indicates  terns  nultiplied  by  sines  and  cosines  of  integral 
multiples  of  £  other  than  those  shown.  The  particular  solutions  corre¬ 
sponding  to  the  terms  shown  on  the  right-hand  sides  of  (4.11)  will  contain 
terms  proportional  to  £  sin  n£  or  £  cos  n£  [n  *  1  for  (4. 11. a,  b)»  n  *  2 
for  (4.11.C,  d)3.  Thus,  the  second  approximation  would  be  unbounded  for 
large  £  while  the  first  approximation  is  bounded  for  all  £.  These 
unbounded  terms  are  called  singular  terms.  The  terms  on  the  right-hand 
sides  of  (4.11)  indicated  by  +  ...  do  not  lead  to  singular  terms. 

The  idea  of  a  perturbation  solution  is  that  higher  order  terms  in 
the  series  solution  represent  small  corrections  to  this  first  term  to 
obtain  a  uniformly  valid  expansion.  The  presence  of  this  singularity 
causes  this  fundamental  idea  to  be  violated.  Therefore,  since  the  expres¬ 
sions  of  n  dependency  are  independent  of  the  variable  causing  the  singu¬ 
larity,  the  n-dependent  expressions  can  be  set  individually  equal  to  zero 
to  avoid  this  problem.  Therefore,  from  equations  (4.11),  the  resulting 
equations,  which  are  eight  ordinary  first-order  differential  equations 
having  n  dependency,  become 

dA 


IT  +  3°A1  +  ^®l  +  ^A1A3  +  B1B3  +  A2A4  +  B2V 


e 


+  ^nidCBjAg  -  AjBg  +  B2A4  -  A^J  *  0 


G 


dB 


T5T  +  +  ’  B1A3  “  A4B2  *  A2B43 


- i  W»  llMIUttfrlf 


♦  “InutA^Ag  ♦  B^Bg  +  AjA4  ♦  BjB43  *  0 

dA*  .  •  « 

HT  *  T°A2  *  #**2  +  2CA1A4  +  B1B4  "  A2A3  ‘  B2B33 
+  |«wC-AxB4  ♦  AUB1  +  A2B3  -  AgBj]  *  0 
dB-  »  % 

TrC  +  rB2  -  rA2  +  ^‘Vl  +  B4A1  -  A2A3  +  A3B23 
+  -iflStA^  ♦  BXBU  -  A3A2  -  B2B33  “  0 

dT  f  2*A3  +  4KB3  *  8  CA22  •  V  +  V  "  V3 

+  ■gwnCA^B^  -  AgBjD  =  0 

dB„  .  , 

JTf2°V  4KA3  +  4CA2B2  "  A131] 

+  ^[A/  -  “  Ax2  +  Bx2]  =  0 

1  1 

g^~  +  2^A4  +  4KB4  +  ^0^2  -  AxA2] 

♦  ■junCA^Bj  ♦  AgB^]  »  0 
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"  ■  -  •  *** w  m 


\r~  '<*4  nnm»^in"'  t  '  -e  j^*i"r.**  T->'”TvrT’ 


•gunCA^j  -  *  0  ♦ 


Since  equations  (4.12)  are  first-order  nonlinear  ordinary  differential 
equations,  the  fourth-order  Runge-Kutta  program,  previously  developed, 
can  be  used  to  solve  for  the  modal  amplitude  coefficients.  By  finding 
these  coefficients  for  various  points  in  time,  a  relation  between  the 
results  of  equation  (3.21)  and  equation  (4.12)  can  be  observed  to  the 
approximation  of  order  e. 

Solving  equations  (4.12)  for  the  highest  derivative  (first  order 
in  this  case)  and  substituting  n  -  et,  the  governing  equations  for  the 
Runge-Kutta  program)  become 


ar  -  ?  (aia3  +  “A  +  A2\  ♦  W 


-  -jnuXBjAg  -  AjB3  +  B2A4  -  A2B4)] 


aBi  i—  1  1 

ST^^B^Pr  ^AIB3  -  B1A3  -  \B2  +  A2B4> 


-  -jnwCA^Ag  +  B^Bg  +  AjA4  +  B2B4)] 


E  C-  'joAj  -  ^^2  -  j(AxA4  +  B^-  AjAg  -  BgBg) 


-  Tgnu)(A4B^  -  AXB4  +  AjBg  -  AgB2)J 


3—  1  1 

3t"  “  e  C“  2°h2  +  2  ^2  "  3*  Vl  "  A4B1  +  A3B2  “ 


a2b3) 


o 


"  Tn(Al\  +  BiB4  -  A3A2  ”  B2B3):1 
dA 

ar  *  «r-  £*3  -  4KB3  -  a<A2!  -  V  ♦  V  -  y> 

-  ■|nu»CA1B1  -  AjBj)] 

dB3  1-  1 

HT*  *  cC~  2°B3  +  4KA3  “  4<A2B2  “  AlV 

-  ^|n«(A22  -  B22  -  Ax2  +  B^)] 

dA  __ 

ar ■  e[-  rA4  -  4KB4  -  l<BiB2  -  AiV 

-  •gna)(A1B2  +  AjB^)] 

dB  _ 

=  e[-  -|oB4  +  4KA4  +  ^32Ai  +  A2B1^ 
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+  fmU^  -  BxB2)] 


(4.13) 


It  is  often  convenient  to  express  the  equations  for  and  in 
terras  of  amplitudes,  C^,  and  phase  angles,  which  are  also  functions 
of  the  slow  time  variable  n.  Mathematically,  wo  can  express  the  relation¬ 
ships  between  the  quantities  as 

Ai  =  Ci  003  ^i  (4. 14. a) 


B^  *  sin 


...) 


(4.14.b) 


tmiliini mi  nl.iiiV, 


L, 


a 


o 
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dA, 

dCi 

i*. 

ZTa 

IT  GOS  *i 

'  °i  3T  ,ln  *1 

(4.\4.c) 

dBi 

dCi 

d^. 

zr  * 

dTj~  sin  ♦i 

+  ci  zr  C08  *i 

(4.14.d) 

where  i  *  1,  2,  3,  and  4  for  each  of  the  equations  above.  Substituting 
the  expressions  of  (4.14)  into  the  first  two  equations  of  (4,12),  the 
resulting  equations  become 


dC 


d$i 


air cos  *i "  ci  zr sln  *1 f  I  cos  *i  + 

+  ^CCj^Cg  cos  ^  cos  $3  *  C^Cj  sin  ^  sin  «p3 

+C2C4  cos  cos  *4  +  C2C4  sin  $2  sin  $4]  +  ^ni^CjCg 

sin  cos  -  Cj.Cg  cos  ^  sin  $3  +  C2C4  cos  <ji4  sin  $2 
-  C2C4  cos  <t>2  sin  $4]  =  0 

dCl  1_  1 

zr sin  *i +  ci  ar  008  *i +  2°ci sin  h  *  f<ci co®  h 

♦  ^ClC3  cos  sin  $3  “  CjCg  cos  sin  <p1 

”  C2C4  008  *4  ®in  $2  +  C2C4  COS  *2  sin  ^  + 

cos  ^  cos  $3  +  C,^  sin  $L  sin  $3  +  C?C4  cos  $2  cos  $4 


O 
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+  C2C4  sin  sin  4>4J  *  0  , 


(4.15) 


Multiplying  th.  first  .guntlon  by  cos  ^  and  tin  second  s,uation  by  sin^, 
adding  the  two  expressions  together,  and  using  appropriate  ■multiple-angle 
identities  from  trigonometry,  the  resulting  equation  for  becomes 

t  "J02!  +  2^Ci<'3^co8^2(^i  ”  ^3^ 

+  C2C4Ccos(«>2  -  *4  +  4»1>3>  +  ^{niC3 

sin(2^  -  d>3>  +  C2C4  sin(#2  -  4>4  +  4^))  *  0  ,  '  VL6) 

Similarly,  multiplying  the  first  equation  of  (4.15)  by  -sin  4^  and  the 
second  equation  by  cos  4^,  adding  the  two  expressions  together,  and  using 
appropriate  multiple-angle  identities  for  trigonometry,  the  resulting 
equation  for  4^  becomes 


d4*  i  ,  ,  C9C„ 

dT  -  ¥  -  ?°3  sln(J\  -  *3>  +  -I7-  sin<W+l)3 

c  c 

+  -2-nwCC3  cos(241-4'3)  +  cos(4»1+4l2-^4)3  =  0  . 


(4.17) 


Using  these  procedures  discussed  above,  equations  for  C2»  $2,  Cg,  4>3>  C4, 
and  4>4  can  be  derived.  Thus,  these  transformed  equations  are 

dC2  1-  1 

3^—  +  ^oC2  +  ^CC^C4  cos(41-4>4+42)  -  C2C3cos(24i2-4i3)] 


+  ^ndj(!CjLC4sin( <t>^-44+4>2 )  -  C2C3sin(24>2-$3)]  =  0 


3T  "  1  K  “  ■in(#1-#4+*a)  -  C3  «in(2*2-*3)J 

c  c 

+  co»(*1-^4+*2)  +  C3ooa(2$2-OJ  *  0 

2 

dC3  1-  1 

HrT  +  2°C3  +  ■£C22°oa(2*2~*3)  “  Vcos^-^)] 

-  l!™CC22sin(2$2-$3)  *  C12sin(2<fr1-^3)]  »  0 

d$3  C  2  c  2 

_  -  4K  +  sin(2^2-<fr3)  -  -i-  sin(2^3)3 

3  3 

C  2  P  2 

1  -  v2  '■'i 

+  ygnwC^—  cos(2$2-*3)  -  cos^-^)]  =  0 

3  3 

dC4  1-1 

dn~  +  2°c4  "  4^C1^2COS  ^ ^  ^  +  in“CClC2 
slnC^  +  -  *4)]  =  0 

1  C1C-)  C.C. 

•**  st-^*  sin(*1+*2-*4):  -  inSt-i-2- 

+  ^2  -  ^4)3  a  0  ,  (4.16 

Equations  (4.16),  (4.17),  and  (4.18)  are  the  general  combustion 
equations  in  terms  of  amplitudes  and  phase  angles.  From  this  point, 
special  cases  can  be  investigated  isolating  certain  conditions  and  closed 
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form  solutions  oan  b«  obtained  for  these  cases.  It  is  convenient  to  do 
this  in  order  to  check  the  closed-form  results  of  the  special  cases  with 
the  results  from  the  general  equations  (4.16),  (4.17),  and  (4.18)  when 
the  same  conditions  are  imposed. 

The  first  case  to  be  evaluated  is  the  case  for  standing  waves 
with  no  combustion  effects.  To  simulate  standing  wave  effect,  set  the 
amplitudes  and  C4  and  phase  angles  ^  and  ^  equal  to  zero.  This 
automatically  satisfies  four  of  the  eight  equations  (4.18).  To  achieve 
the  no-combustion  effect,  set  the  interaction  index,  n,  equal  to  zero. 
Also,  set  the  correction  variable,  K,  equal  to  zero  since  the  effect  of 
K  will  be  investigated  separately  at  a  later  time.  Imposing  these  con¬ 
ditions,  the  governing  equations  reduce  to 

dC .  i  i 

oC^  +  cos(2$;L-^3)  =  0  (4. 19. a) 

d$i  i 

J  [ CgSin^-.frg)]  =  0  (4.19.b) 


dC3  1  -  1 

3T  +  2  oC3  ’  8  C12cos(2^1-4.3)  =  0 
d$3  1  Cl2 

3 GT  "  8C7~  ain(2V*3)  *  0# 

The  initial  conditions  imposed  for  this  case  are 

cx(o)  =  1 


(4.19.c) 

(4.19.d) 


u 


C3(0)  *  0 


♦l<°>  ■  *10 


♦3(0>  *  ♦so  . 


(4.20) 


To  attempt  a  closed-form  solution,  lat 


C1  *  ®’^r'Fi 


(4. 21. a) 


C3  *  •'%5nr3 


(4.21.b) 


-  dF, 


♦  •"*V) 


(4.21.C) 


dC„  , -  ,  -  dF- 

3^  .  5>r3  ♦  .-*V»  . 


(4.21.d) 


Substituting  these  expressions  into  equations  (4. 19. a)  and  (4.19.c)  ana 
dividing  through  by  e"^an,  the  resulting  equations  become 


i  t  \  cos(2^  -  *3)e-*5riFlF3  =  ° 


(4. 22. a) 


~  -  •§•  cos^  -  ♦3)e"%5nF12  =  0  . 


(4.,  22. b) 


Multiplying  equation  (4.2  2.a)  by  1/4  and  equation  (4.2  2.b)  by  F^/F^^  and 

wLiin 

adding  the  two  equations,  terras  containing  the  cos(2<^  -  $j)e  are 
eliminated.  In  doing  so,  the  result  becomes 


Q 


Letting  2$^  -  $3  =  £tr,  which  satisfies  equations  (4.19.b,  d),  yields 
003(24*2  -  4>3)  *  (-1)  where  l  *  0,  1,2,3.  .  .  Substituting  this  expres¬ 
sion  and  integrating  the  above  equation,  the  resulting  equation  becomes 


— !7H?VIT>  ,'T’  '',(«WCT!",,,*'^T''TT  ’I"  !,*■  ' '  V'  ”  ■ 


\  tanh*1  2Fg  *  C-  |  e^+D^C-l)* 


(4.28) 


where  is  a  constant  of  integration.  Using  the  initial  condition  F3(0)» 
0,  then,  it  can  be  shown  that  D2  *  2/a.  Substituting  and  taking  the 
hyperbolic  tangent  of  both  sides  of  equation  (4.28),  the  result  becomes 


F3  *ljtanhC|-<-lAl-e-%5n)3  , 


(4.29) 


Substituting  this  expression  into  equation  (4.26)  and  simplifying,  the 
resulting  equation  becomes 


r  .  sech[^  (l-e“*°n)3  . 
x  2o 


(4.30) 


Substituting  equations  (4.29)  and  (4.30)  into  equations  (4. 21. a)  and 
(4.21.b),  and  substituting  n  =  et  and  u>  =  5e,  the  resulting  closed-form 
solution  for  wave  amplitudes  and  C3  are 


C  *  e^fsechC^i^l-e^)]} 

A  2u 


(4. 31, a) 


i  C 


C3  a  { tanhC^- l-e~^it ) 3  )  , 

3  2  2u 


(4.31.b) 


To  find  expressions  for  and  $3,  substitute  the  relation  that  2$,-$3s£it 
into  equations  (4.19.b)  and  (4.19.d)  and  integrate  and  evaluate  the  con¬ 
stants  of  integration  with  the  initial  conditions;  the  results  are 


*1  =  *10 


I 


IT  WWIWW  "V»p v*r,  ' T 


p;  «n«r»*  jwr*  rwr*'*  t  wp*** 1 


h*  ho*  2 ho 
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(.4.32) 


where  is  a  constant  and  $>3  is  bn  radians  out  of  phase  with  2^,  It 
can  be  seen  that  a  special  set  of  initial  conditions  is  necessary  to  be 
consistent  with  this  solution.  A  representative  set  is  <^0  *  $gQ  *  0 
which  corresponds  to  L  *  0. 

Inspection  of  equations  (4.31)  reveals  that  the  magnitude  of 
continually  decreases  with  time  while  the  magnitude  of  Cg  first  increases 
and  then  decreases.  An  interesting  special  case  of  equations  (4.31) 
occurs  in  the  absence  of  steady-otate  combustion  (»  =  0) .  The  results  of 
this  case  are 


Cx  =  sech[--~-'~— ■"] 


C3  =  j  tanh[(-fe- ] 


A 


b 


(4.33) 


These  results  show  that  a  disturbance  in  the  form  of  the  first  mode  is 
transferred  to  the  second  mode  as  time  increases.  It  is  thought  that  this 
indicates  the  beginning  of  the  steepening  that  leads  to  the  formation  of  a 
shock  wave.  It  can  be  seen  that  the  presence  of  damping,  in  the  form  of 
steady-state  combustion,  inhibits  this  process. 

The  second  case  to  be  investigated  is  that  of  standing  waver  with 
gas-dynamic  nonlinearities  neglected.  To  simulate  the  standing  wave 
effect,  let  the  amplitudes  C2  and  and  the  phase  angles  <J>2  and  equal 
zero.  Again,  this  automatically  satisfies  four  of  the  eight  equations  of 
(4.18).  To  achieve  omission  of  gas-dynamic  nonlinearities,  let  i  =  0. 
Also,  let  the  correction  variable,  K,  be  equal  to  zero  for  simplicity. 
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In  doing  so,  the  rssulting  aquations,  based  upon  equation  C4.18), 
become 


1  1  -  i  _ 

3n“  +  2“  +  1  niJ^ciC3  8^(2^  -  $3)J  8  0 

^1  i 

dT+2  nuCC3  c°s(2*i  -  *3)H  s  0 
dC3  1  1 

dT  +  2  ^3  "  16  sin<2*l  -  *3)]  *  0 

d*3  I  -  Cl2 

dlT  +  16  ntoC“c"  Cos(2*l  “  *3)]  =  0  •  (4*34) 

3 

The  initial  conditions  imposed  for  this  case  are 
(^(O)  a  1 

c3(o)  =  0 

*1<°>  *  ho 

♦3(°)  =  $30  .  (4.35) 

ket  2<^1  ”  ^3  a  ^2^  +  ^  =  0*  i»  2  .  .  .  .  This  implies  that  sin^^ 

-  $3)  =  (-1)  and  cosO^  -  $3)  =  0.  Substituting  into  (4,34)  and  solving 
in  the  manner  indicated  previously  one  obtains  expressions  for  the  ampli¬ 
tudes  for  and  C3  which  are 

Cx  =  e'^tsecC^  net-lAl-e"^1)]} 


(4. 36. a) 


I  **V"«j»r,<'i*i  w«^*iw  vw'A'»-T'’¥*»r*  "tm* '  ■v-i’ 


^-%nt£  ry^  i\^/i  »%wt% 

C3  a  m*2J2  *  4  R&(jmlj  CL*»0  ) 


*1 a  ♦« 


♦3  *  J*io  -  <iTi>* 


(4,36.b) 


(4.36.c) 


(4,36.d) 


-her.  1.  constant  and  ♦,  i,  (2tnM  radlans  out  of  phase  ^  * 

A*  In  the  previous  solution,  speciai  initial  conditions  are  required  to 
produce  this  solution.  A  representative  set  is  *1Q  *  0,  ♦3Q  *  -*/2,  which 
corresponds  to  l  =  0. 

The  secant  and  tangent  both  become  infinite  when  their  arguments 
take  on  the  value  *»/2.  In  (4.36. a,  b),  the  arguments  of  these  functions 
start  at  zero  at  t  *  0  and  have  a  maximum  absolute  value  at  ne/23/2. 

,  if  nc/2  <  s/2,  the  tangent  and  secant  never  become  infinite  and 

Cl  and  °3  6''<”’tUaUy  d*Cay  to  “«  4“  *>  «■.  influence  of  the  exponential 

function.  This  is  a  stable  situation.  If.  on  the  other  hand,  ns/23'*  > 

*/2,  the  tangent  and  secant  become  infinite  at  t.  .  (2/u>|f  ntl-2%v/(nt):i| 

causing  Cx  and  C3  to  become  infinite.  This  is  an  unstable  situation. 

Thus,  the  boundary  between  stable  and  unstable  behavior  is  indicated  by 
the  equation 


ne/2^^  *  ir/2  . 

The  stability  equation  in  the  n-e  plane  has  the  form 


(4.37) 


n  *  2Vc  =  4.442/e  , 


(4.38) 


This  has  the  form  of 


a  rectangular  hyperbola  and  is  independent  of  Z. 
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For  the  case  of  traveling  waves ,  it  is  more  convenient  to  work 
with  the  general  perturbation  equations  expressed  in  modal  amplitudes 
in  terms  of  the  real  time  variables,  equation  (4.13).  To  simulate  the 
effect  of  spinning  or  traveling  waves,  let  the  following  modal  ampli¬ 
tudes  be  equal.  These  relations  are 


B 


4 


A 


3 


B3  “  “A4  .  (4.39) 

It  can  be  shown  that  substituting  the  relations  (4.39)  into  equation 
(4.10),  expressing  the  results  in  terms  of  the  real  time  variables,  sub¬ 
stituting  these  expressions  into  equation  (3.18),  and  using  appropriate 
multiple-angle  formulas  leads  to 

♦(0,t)  s  A1cos(t-0)  -  A2sin(t-9)  +  AgCOS  2(t-0) 

-A^sin  2(t-0)  +  .  .  .  .  (4.40) 

which  has  the  form  of  a  sura  of  traveling  waves.  Substituting  the  expres¬ 
sions  in  (4.39)  into  equations  (4.13),  these  eight  equations  reduce  to 
four  pairs  of  identical  equations.  The  four  independent  equations  listed 
below  are 


^A1  1-1 

-  eC-  oA1t^KA2-i(A1A3+A2A4)-nw(A1A4-A2A3)] 


**«.  .1  *« **«*3»*  «*••»?*•  ■«■*«*  Nt  ■*  > alX*S*W«v'r  >«v*.w 


1-1 

“  S  eC“  2  oA2  "  2  KAri(AlA4"A2A3)+nu,CA2A4+A1A3)] 


2  4  13 


et-  oA3+4KA4+  ^i(Aj^-A2a)+^nw(A.jA2)3 


dAu  .  .  . 

*  cC-  7  5A4-4KA3+  i  KA^j)-  |nS(A12-A22)]  . 


(4.41) 


By  making  the  substitution,  we  have  reduced  to  a  system  of  four  equations 
and  four  unknowns.  By  solving  for  the  modal  amplitudes  A^,  the  modal 
amplitudes  are  readily  computed  by  using  the  relations  of  (4.39)  to 
determine  the  entire  nature  of  the  wave  form. 

For  the  case  of  traveling  waves  omitting  gas-dynamic  nonlinearities, 
let  the  amplitudes  A^  and  Ag  equal  zero.  Then  set  i,  the  gas-dynamic 
index,  equal  to  zero.  Again,  for  simplicity,  let  the  correction  variable, 
K,  controlling  physical  chamber  configurations,  be  zero.  In  doing  so,  in 
terms  of  the  transformation  variable,  n»  the  resulting  equations  become 


w2  1  -  - 

dir*  +  2  °  A2  nwCAjA^]  =  0 


uo4  1  -  1  - 

3T  +  ?#a»  -7n«A2  ■  “ 


(4,42) 


which  is  a  system  of  two  equations  and  two  unknown  modal  amplitudes.  To 
find  an  exact  closed-form  solution  to  these  equations,  let 


a  =  “%orip 
a2  -  e 


\  •  «-i5nF. 
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Fx  *  8ecc44(l-‘^n)3 


„  _  tarjr*^  n«  ,,  _-%orKi 

F0  »  —  (l-«  )J  . 

2  2^2  Z  5 


(4.46) 


Expressing  the  results  of  (4.45)  in  terms  of  modal  amplitudes  by  substi¬ 
tuting  into  (4.42),  the  resulting  equations  become 


a.  .  i§ 

4  4  a 


k-%on 


A4  a 


2 - tant—2^  d-e^00)] 

2/2  a 


(4.47) 


The  results  for  traveling  waves  (4.47)  are  quite  similar  to  the  results 
for  standing  waves  (4.36)  for  the  case  of  no  gas-dynamic  nonlinearities. 
The  same  behavior  can  be  expected  as  was  dircussed  in  the  standing  wave 
cane  about  the  nature  of  oscillation  of  the  modal  amplitudes.  The  only 
significant  difference  is  the  value  to  determine  the  boundary  of  stability 
for  the  interaction  index  governing  the  combustion  terms.  The  stability 
condition  for  traveling  waves  is 


/2  _  * 
"2  n  2 


Thus,  the  equation  of  the  stability  boundary  in  the  n-e  plane  is 


w  2.22 

n  -  — tt*  =  "■ 

2*e  e 


(4.48) 


(4.49) 


Comparing  equation  (4.49)  to  (4.38)  shows  that  the  stability  boundary  for 
the  interaction  index  is  half  as  great  for  the  traveling  wave  case  as  for 
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the  standing  wava  case  for  any  c.  This  will  be  verified  in  a  later  pre¬ 
sentation  of  results  of  various  numerical  cases. 

For  the  case  of  traveling  waves  with  no  combustion,  let  the  ampli¬ 
tudes  A2  and  A4  equal  to  aero.  Then  set  n,  the  interaction  index,  equal 
to  aero,  and,  again,  let  the  correction  variable  K  equal  to  zero.  Sub¬ 
stituting  into  equations  (4.40)  and  transforming  into  variable  n*  the 
results  are 


u/vl  i  . 

dT+2aV  A1A3  3  ° 


dA3  1  -  A  1  2 

ItT  ?  0  A3  “  4  A1 


(4.50) 


with  initial  conditions 


A1(0)  =  1 


a3(o)  =  0 


which  again  is  a  system  of  two  equations  and  two  unknown  modal  amplitudes. 
To  find  an  exact  closed-form  solution  to  these  equations,  use  similar 
transformations  as  shewn  in  (".42).  In  doing  so,  and  simplifying,  the 
results  are 


1  +  a~b°ri  _  _ 

5T  *  8  V2  '  0 


i .  A  a  ,  o 


m»  «***• 

dn  4 


(4.51) 
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with  initial  condition* 


F^O)  *  1 


Fj(0)  *  0  , 


Solving  those  equations  in  the  same  manner  as  before,  the  results  are 


i. 


Fa  *  sech[l/5(l-e"%an)3 
F2  =  i  tanh[l/5(l-e"^n)]  . 


(4.52) 


Again,  expressing  the  results  of  (4.51)  in  terms  of  the  modal  amplitudes 
of  the  form  of  equation  (4.43),  the  resulting  equations  become 


A^  *  e~^°nsechCl/5(l-e”^n)] 


~%cm  i  ■* 

A3  «  2_  tanh[l/o(  l-e~*on)D  . 


(4.53) 


The  results  for  the  traveling  waves  (4.52)  are  similar  to  the  results 
for  standing  waves  (4.31)  for  the  case  of  no  combustion.  A  disturbance 
initially  having  the  form  of  the  first  mode  eventually  is  transformed  into 
one  having  the  form  of  the  second  mode.  To  compare  these  results  for 
standing  waves  and  traveling  waves  to  the  general  perturbation  equations, 
two  computer  programs  were  written  (Appendices  D  and  E)  which  numerically 
evaluate  the  modal  amplitudes  of  various  conditions  for  standing  and 
traveling  waves. 
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One  last  special  case  is  an  investigation  of  the  effect  of  the 
correction  variable  K.  In  the  special  cases  previously  discussed,  the 
correction  variable  K  was  set  equal  to  zero.  But,  in  this  discussion, 
the  correction  variable  K  will  be  of  primary  importance  in  the  equations. 
To  start  this  analysis,  refer  to  equations  (3.21).  Based  upon  these 
equations,  impose  the  following  conditions.  First,  neglect  combustion 
effects  (i.e.,  n  *  0).  Then,  let  us  consider  only  the  case  of  standing 
waves  (i.e.,  g^  *  g ^  =  0).  Finally,  let  us  neglect  the  steady  state 
burning  rate  (i.e.,  5=0)  and  assume  that  the  terms  multiplied  by  eK 
are  larger  than  those  multiplied  by  e  above.  This  can  be  accomplished 
by  writing 

Kx  =  eK  (4.54) 

and  treating  as  a  quantity  of  0(1).  Imposing  the  above  conditions 
and  substituting  equation  (4.54)  into  the  equations  (3.21),  the  result¬ 
ing  equations  become 


d2f. 

df  df 

[1+y 

i 

dt2 

*  fl  +  2eCf2  lt~  +  fl  dt""-1  *  0 

(4. 55. a) 

d2f2 

dfl 

[1+4^] 

£ 

dt2 

'  +  4f2  -  efl  dT  "  0 

(4.55.b) 

with  initial  conditions 


fi'°>  *  1 
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dfl 

^-CO)  *  0 


f2  <o>  *  0 

cft^O)  *  0  . 


First i  assume  a  straightforward  perturbation  solution  similar  to  the 
equations  (4.6)  except  the  functions  are  dependent  upon  the  real  time  t. 
Substituting  these  assumed  solutions  into  the  equations  and  initial  con¬ 
ditions  of  (4.55)  and  keeping  terms  of  0(1)  and  0(e),  the  separated 
equations  become 


d2f 


10 


dt2 


■  JL 

+  Ti+k^T 


k10 


=  0 


(4. 56. a) 


d2f, 


20 


dt2 


(1+4KX) 


20 


*  0 


(4.56.b) 


d2f 


11 


dt2 


+  ( 


■)f. 


1+K, '11 


Trap1-* 


df10  .  df20l 
2o“dt“  ‘  holT1 


(4.56.c) 


d2f21  4 

_ £L  +  (_ _ >f 

4t*  lt4Kl  21 


nwk"yCfi<rara 


( 4. 56.d) 


with  initial  conditions 


f10«»  ■  i 


fu(0)  »  0 
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df 

-—-—(0)  *  0 


dfn 
— jynO) 


»  0 


f20(o)  *  0 


(0)  =  0 


df20(0) 

dt 


a  0 


■  df*  .(0) 

xx -  a  0  , 


The  first-order  equations  (4. 56. a  and  b)  can  be  solved  by  assuming  the 
usual  assumed  solution  for  linear  differential  equations.  Doing  this 
and  applying  the  appropriate  initial  conditions,  the  results  for  the 
first-order  terms  are 

fio  ■  003JiS^  1 


(4.57) 


Substituting  (4.57)  into  the  right-hand  side  of  (4.56.c)  the  equation 
becomes  a  homogeneous  linear  differential  equation.  Solving  in  the 
usual  manner  and  applying  the  appropriate  initial  conditions 

f1L  =  0  .  (4.58) 

Substituting  (4.57)  into  the  right-hand  side  of  equation  (4.56.d),  the 
resulting  equation  becomes  a  linear  differential  equation  with  a  particu¬ 
lar  solution.  By  assuming  an  appropriate  homogeneous  and  particular 
solution  and  evaluating  the  constants  using  the  appropriate  initial  condi¬ 
tions,  the  result  becomes 


<4.59) 


must  be 
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A1<n) 

Bl(n) 

By  similar  procedure,  it  can 
and  comparing  it  to  equation 

A3(n> 

V> 

To  show  the  validity  of  equations  (4.63)  and  (4.64),  the  problem  is  now 
solved  using  equations  (4.12)  which  are  derived  from  equations  (3.21)  by 
the  use  of  the  two- variable  perturbation  method.  To  reproduce  the  condi¬ 
tions  imposed  on  the  problem  just  discussed,  let  there  be  no  combustion 
(i.e.,  n  =  0),  let  there  be  no  steady-state  burning  rate  (i.e.,  a  =  0), 
and  let  there  be  only  standing  waves  existing  (i.e.,  A2  =  *  B2  *  * 

0),  Imposing  these  conditions  on  equations  (4.12),  the  resulting  equa¬ 
tions  become 


a  cos  j  Kn 


=  sin  j  Kn  . 


(4.63) 


be  shown  that  evaluating  equation  (4.60.b) 
(4,61.b),  the  results  are 


=  2^Csin  4Kn  -  sin  Kn] 


24K1 


;cos  Kn  -  cos  4Kn]  . 


(4.64) 


nr  ~ 


dn 


+ 


4KB3 


+  ■jCA^Ag  +  B^iJg]  = 

+  |tAlB3  *  BlV  * 
+  fci*  -  * 0 


0 


0 


.i L 
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than  substituting  these  expressions  into  the  equations  C4.66)  and  keeping 
terms  of  0(1)  and  0(1/K),  the  resulting  separated  equations  become 


dA,  _  . 

~dF"  +  2  Bio  8  0 


(4. 68. a) 


dB1n  . 

iQ  _  i  a  s  0 

TET  2  A10  0 


(4.68.b) 


^30 

TT  +  4B3o  s  0 


(4.68.c) 


^B30 

■dT  ■  4A30  =  0 


(4.68.d) 


|  i 

~  2  n 


“  2^A10A30  +  B10B30-* 


C4.6S.e) 


dB 


— ii  -  i  a 
dc  2  11 


“  2tA10B30  '  B10A30^ 


(4.68.f) 


dA 


d?  +  4B31  *  “  8^B102  “  A102;3 


(4.68.g) 


dB 


d?  '  4A31  *  4^A10B10^ 


(4.68.h) 


with  the  initial  conditions 


V°>  ■  1  B10<0>  *  0 


AU<°)  *  0 


Bn(0)  =  0 
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A30(0)  “  0  B3Q(0)  "  ° 

A31(0)  =  0  B31(o>  a  0  * 


Since  the  first-order  equations  are  coupled,  differentiate  equations 
(4. 68. a  and  c)  once  with  respect  to  c  then  substitute  equations  (4.68,b 
and  d)  into  these  equations  resulting  in 


d*A 


d?2 


+  1  A 
4  A10 


0 


d?2 


+  16A 


30 


(4.69) 


As  can  be  seen,  equations  (4.69)  are  linear  differential  equations 
which  can  be  evaluated  by  the  usual  manner.  In  doing  so  and  applying 
the  appropriate  initial  conditions,  the  resulting  first-order  modal  ampli¬ 
tudes  are 


A1q  »  cos  *  cos  -jKn 

A30  .  0  .  (4.70) 

Knowing  values  for  A1Q  and  A3Q,  substitute  these  values  into  equations 
(4.68.b  and  d)  and  apply  appropriate  initial  conditions.  The  results 
become 

B1q  *  sin  =  sin  -jKn 


(4.71) 


•  v**r  ' 
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Substituting  the  results  of  (4.70)  and  (4.71)  into  the  right-hand  side  of 
equations  (4.68.e-h),  the  resulting  equations  become 


dAn  . 
d? 


T  Bii 


^11  .  1 
-  - A 


2  °11 


(4. 72. a) 


(4.72.b) 


+  4B 


31  8 


cos  c 


(4.72.C) 


4A, 


31 


sin  C  , 


(4.72.d) 


Since  equations  (4.72.C  and  d)  are  coupled,  differentiate  both  equations 
once  with  respect  to  c  and  substituting  equations  (4.72.C  and  d)  into  the 
appropriate  terms  of  the  new  set  of  equations,  the  resulting  equations  are 


d^A 


31 


dC2 


+  16A 


31 


sin  c 


d*B 

dC 


31 

2 


+  16B31 


5 

8  C0S 


C 


(4.73) 


Equations  (4,73)  are  a  set  of  linear  differential  equations  with  homo¬ 
geneous  and  particular  solutions.  Solving  these  equations  in  the  usual 
manner  and  using  the  appropriate  initial  conditions,  the  resulting  modal 
amplitudes  are 


A31  8  ■j—  (sin4e  -  sine)  =  -^(sin4Kn  -  sinKn) 


, . *»*1  JgRW!  ii  «»*k, 
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B31  *  2^cos<*  “  cos4c)  =  ~<cosKn  -  oos4Kn).  (.4.74) 

In  a  similar  manner,  the  results  for  the  modal  ampHxudes  and  can 
be  determined  to  be 


B11  =  0  <4-75' 

evaluated  with  the  appropriate  initial  conditions.  Therefore,  substitut¬ 
ing  the  results  of  (4.70),  (4.71),  (4.73),  and  (4.74)  into  the  assumed 
perturbation  solution  of  (4.67),  the  resulting  modal  amplitudes  become 

=  cos  4  Kn  +  .  .  . 

=  sin  y  Kn  +  .  .  . 

Ag  *  ^~^<sin  4Kn  -  sin  Kn)  +  .  .  . 

Bg  =  ~£<3in  4Kn  -  sin  Kn)  +  .  .  .  .  (4.76) 

It  can  be  seen  that  equations  (4.76)  are  identical  to  equations  (4.63)  and 
(4.64).  This  indicates  that  the  two-variable  method  produces  the  correct 
solution.  Equations  (4.60)  indicate  that  the  presence  of  K  changes  the 
frequency  of  each  of  the  first  two  acoustic  modes  and  further  renders  the 
ratio  of  the  second  frequency  to  the  first  a  non-integer  number  in  general. 
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Equations  (4.76)  show  how  this  offset  manifests  itself  in  the  two-variable 
perturbation  solution. 

These  results  can  be  used  in  another  way.  If  the  nonlinear  terms 
are  neglected  in  (4. 5 5. a),  the  results  are 


(1+Kp 


+  f. 


3 


0 


d2f„  df. 

(I+4K  )  - i  +  4f.  -  ef.  ~  =  0 

1  dt2  2  1  dt 


fx(0)  =  1, 


dfx(0) 

dt 


=  0,  f2(C)  =  0, 


df2(0) 

dt 


=  0 .  (4.77) 


It  c&  y  shown  tha.  equations  (4.60)  constitute  the  exact  solution 
of  equation  (4.77),  If  the  corresponding  terms  are  neglected  in  equations 
(4.65),  the  results  are 


0 


0 


dA3  1 

ir +  4KB3  +  8<B12  “  Ai2)  *  0 


dB,  , 

—  -  4KA.  -  ~  A. B.  *  0 
dn  3  4  11 


(4.78) 


where 


74 


Aj.tQ)  a  1 
B1(0)  =  0 

a3(o)  *  0 
b3(o)  =  0  , 

It  can  be  shown  that  equations  (4.76)  are  the  exact  solution  of  equation 
(4.78).  These  facts  ware  used  to  check  the  accuracy  of  the  computer  pro¬ 
grams  to  be  discussed  later. 

In  the  remainder  of  this  thesis*  a  comparison  of  the  magnitudes 
of  the  modal  amplitudes  will  be  represented  in  graphical  and  tabular 
form.  Under  a  given  set  of  conditions,  the  acoustic  modal  amplitude  pro¬ 
gram,  the  general  perturbation  program,  and  the  analytical  cases  that 
were  programmed  will  be  used  and  results  compared.  Varying  certain  con¬ 
ditions  will  show  their  effect  on  the  changes  in  magnitude  of  the  modal 
amplitudes  through  a  set  range  of  time  which  is  related  to  maintaining 
stability.  From  these  various  cases,  it  will  be  determined  which  param¬ 
eters  and  conditions  have  the  greatest  effect  in  changing  «w.hl  ampli¬ 
tudes  and  which  in  turn  affect  the  stability  criteria  for  combustion 
by  the  methods  discussed  above. 


Chapter  5 


DISCUSSION  AND  PRESENTATION  OF  RESULTS 

In  this  chapter,  results  are  presented  both  in  graphical  and 
tabular  form  which  are  representative  of  the  results  generated  by  the 
programs  listed  in  the  Appendices  B  through  E.  From  these  representative 
sets  of  results,  basic  observations  will  be  made  to  observe  which 
parameters  or  conditions  have  the  greatest  effects  on  the  problems  of 
stability. 

In  Figures  3  and  4,  modal  amplitudes  F^  and  F2  are  graphically 
represented  versus  time  for  a  stable  standing  wave  case.  For  these 
figures,  F1(0)  =  0,  F1’(0)  =  1,  F2<0)  *  0,  F^CO)  =  0,  G1(0)  =  0,  G1’(0)  = 
0,  G2(0)  =  0,  G2'(0)  =  0,  n  =  35,  i  =  1,  K  =  0,  e  =  0.1  and  i  =  0.1. 

The  step  size  used  was  0.1.  Experimentation  showed  that  this  was  a  small 
enough  step  size  to  produce  accurate  results  and  was  used  throughout. 

From  these  figures,  one  notices  that  both  the  first  and  second  order  modal 
amplitudes  decrease  in  amplitude  with  increasing  time.  Also,  F2,  the 
second  order  modal  amplitude,  tends  to  oscillate  at  twice  the  frequency  of 
F^.  These  figures  are  based  upon  one  set  of  parametric  values;  however, 
these  figures  represent  qualitatively  the  results  obtained  using  a  wide 
variety  of  initial  conditions  and  parametric  values.  In  Figures  5  through 
8,  modal  amplitudes  Fj,  F2,  Gj,  and  G2  are  graphically  represented  versus 
time  for  a  stable  traveling  wave  case.  For  these  figures,  F^(0)  =  0, 

Fj*  (0)  ■  -1,  F2(0)  *  0,  G^O)  *  1,  G1 '  (0)  =  0,  G2(0)  =  0,  G2’(0)  =  0, 
n  *  15,  1  «  1,  K  =  0,  u  a  0.1  and  e  =  0.1.  The  general  shape  of  the 


curves  and  the  relative  frequencies  of  oscillation  are  qualitatively 
similar  to  the  stable  standing  wave  case. 

In  Figures  9  and  10,  modal  amplitudes  and  F2  are  graphically 
represented  versus  time  for  an  unstable  standing  wave  case  with  the 
same  conditions  as  the  stable  case  except  that  n  s  50,  As  can  be  seen, 
the  maximum  amplitude  of  F^  starts  to  decrease  then  increase  dramatically 
for  increasing  time.  The  maximum  amplitude  of  Fj  increases  continuously. 

Ih  Figures  11  through  14,  modal  amplitudes  F^,  Fj,  6^,  and  Gj  are 
represented  versus  time  for  an  unstable  traveling  wave  case.  Again,  the 
conditions  are  the  same  as  for  the  stable  traveling  wave  case  except  that 
n  =  30.  Drastic  increases  in  amplitudes  are  observed  for  all  the  modal 
amplitudes  shown  as  time  increases.  The  behavior  is  similar  to  the 
unstable  standing  wave  case.  The  period  of  time  for  traveling  waves  to 
become  unstable  is  about  one-half  the  period  of  time  for  standing  waves 
to  become  unstable.  Thus,  it  seems  that  traveling  waves  are  less 
stable  than  are  standing  waves. 

In  Tables  1  and  2,  a  comparison  of  results  is  presented  for  modal 
amplitudes  F^  and  F2  for  a  stable  standing  wave  case.  For  these  cases, 
F^Q)  =  0,  Ft'  (0)  =  1,  F2(0)  =  0,  F2’(0)  =  0,  G^O)  =  0,  G1'(0)  =  0, 

G?(0)  =  0,  G2’(0)  =  0,  n  =  60,  e  =  0.1,  and  u>  *  0.1.  These  tables 
quantitatively  show  the  effect  of  neglecting  gas  dynamic  non-linearities 
on  the  accuracy  of  the  computations.  Also,  a  comparison  can  be  made 
between  the  exact  solution  method  (Appendix  B  program)  and  the  perturbation 
solution  method  (Appendix  C  program).  From  Table  1,  one  can  observe  that 
the  effect  of  neglecting  gas-dynamic  nonlinearities  is  small  where 
quantitatively  comparing  values  of  the  modal  amplitude  Even  though, 

quantitatively,  the  values  for  the  exact  solutions  and  perturbation 
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Tabl«  1.  Comparison  of  Results  for  Fx  8howing  Effects  of  G«s  Dynamic 
Index  (i)  -  (Fjl  •  0,  Fj,'  ■  1,  Fj  »  0,  F2*  *  0,  »  0,  G1’  ■  0,  G2  *  0, 

G2'  «  0)  -  Stable  Cases  (n  ■  60)  «  Standing  Haves 


i  ■  i 
K  *  1 


Exact 

Solution 

0.19699 

0.36335 

0.55252 

0.69856 

0.81627 

0.90132 

0.95043 

0.96159 

0.93432 

0.86985 

0.77125 

0.64330 

0.49211 

0.32469 

0.14827 

-0.03017 

-0.20430 

-0.36859 

-0.51829 

-0.64917 

-0.75738 

-0.83921 

-0.89118 

-0.91029 

-0.89444 

-0.84300 

-0.75726 

-0.64071 

-0.49885 

-0.33868 

-0.16794 

0.00574 

0.17554 

0.33582 

0.48221 

0.61134 

0.72042 

0.80682 

0.86778 

0.90042 


Perturbation 

Solution 

0.18712 

0.36426 

0.52540 

0.66518 

0.77905 

0.86340 

0.91572 

0.93461 

0.91986 

0.87244 

0.79443 

0.68895 

0.56003 

0.41247 

0.25167 

0.08340 

-0.08633 

-0.25158 

-0.40659 

-0.54604 

-0.66519 

-0.76006 

-0.82756 

-0.86556 

-0.87297 

-0.84979 

-0.79704 

-0.71679 

-0.61200 

-0.48646 

-0.34466 

-0.19159 

-0.03259 

0.12685 

0.28126 

0.42538 

0.55435 

0.66387 

0.75032 

0.81090 


Exact 

Solution 

0.19699 

0.38336 

0.55259 

0.69885 

0.81719 

0.90354 

0.95489 

0.96936 

0.94632 

0.88656 

0.79242 

0.66783 

0.51822 

0.35021 

0.17115 

-0.01142 

-0.19032 

-0.35912 

-0.51242 

-0.64583 

-0.75583 

-0.83953 

-0.89450 

-0.91868 

-0.91049 

-0.86909 

-0.79483 

-0.68959 

-0.55708 

-0.40279 

-0.23365 

-0.05742 

0.11809 

0.28578 

0.43974 

0.57552 

0.68999 

0.78099 

0.84693 

0.88644 


1*0 
K  *  1 

Perturbation 

Solution 

0.18702 

0.36386 

0.52462 

0.66400 

0.77758 

0.86186 

0.91443 

0.93399 

0.92040 

0.87466 

0.79884 

0.69602 

0.57016 

0.42593 

0.26856 

0.10366 

-0.06300 

-0.22566 

-0.37879 

-0.51727 

-0.63655 

-0.73278 

-0.80297 

-0.84504 

-0.85789 

-0.84143 

-0.79656 

-0.72514 

-0.62988 

-0.51427 

-0.38244 

-0.23901 

-0.08892 

0.06271 

0.21080 

0.35042 

0.47703 

0.58653 

0.67550 

0.74121 
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Table  2.  Comparison  of  Results  for  F2  Showing  Effects  of  Gas  Dynamic 
Index  (i)  -  (Fx  ■  0,  Fj'  *  1.  F2  ■  0,  F2*  «  0,  Gj,  ■  0,  G1'  »  0,  Gj  *  0, 
Gj'  »  0)  -  Stable  Cases  (n  *  60)  -  Standing  Waves 


7Z  fr¬ 
it  11 

i  »  0 

K  *  1 

t 

Exact 

Perturbation 

Exact 

Perturbation 

Solution 

Solution 

Solution 

Solution 

0.2 

0.00012 

-0.00485 

0.00003 

-0.00253 

0.4 

0.00113 

-0.01309 

0.00043 

-0.00938 

0.6 

0.00422 

-0.02223 

0.00205 

-0.01865 

0.8 

0.01060 

-0.02944 

0.00602 

-0.02784 

1.0 

0.02110 

-0.03215 

0.01336 

-0.03424 

1.2 

0.03582 

-0.02854 

0.02471 

-0.03553 

1.4 

0.05397 

-0.01795 

0.03497 

-0.03023 

1.6 

0.07375 

-0.00104 

0.05816 

-0.01798 

1.8 

0.09260 

0.02022 

0.07742 

0.00026 

2.0 

0.10749 

0.04283 

0.09521 

0.02232 

2.2 

0.11543 

0.06317 

0.10863 

0.04514 

2.4 

0.11407 

0.07764 

0.11493 

0.06516 

2.6 

0.10215 

0.08319 

0.11198 

0.07891 

2.8 

0.07988 

0.07788 

0.09878 

0.08357 

3.0 

0.04909 

0.06128 

0.07573 

0.07745 

3.2 

0.01309 

0,03461 

0.04478 

0.06031 

3.4 

-0.02375 

0.00071 

0.00929 

0.03352 

3.6 

-0.05653 

-0.03632 

-0.02639 

-0.00006 

3.8 

-0.08051 

-0.07164 

-0.05749 

-0.03640 

4.0 

-0.09180 

-0.10031 

-0.07945 

-0.07083 

4.2 

-0.08798 

-0,11801 

-0.08865 

-0.09863 

4.4 

-0.06850 

-0,12165 

-0.08296 

-0.11571 

4.6 

-0.03490 

-0.10993 

-0.06211 

-0.11921 

4.8 

0.00924 

-0.08353 

-0.02786 

-0.10793 

5.0 

0.05868 

-0.04515 

0.01610 

-0.08257 

5.2 

0.10711 

0.00077 

0.06458 

-0.04571 

5.4 

0.14796 

0.04864 

0.11144 

-0.00152 

5.6 

0.17534 

0.09236 

0.15041 

0.04468 

5.8 

0.18493 

0,12615 

0.17593 

0.08715 

6.0 

0.17466 

0.14534 

0.18393 

0.12042 

6.2 

0.14513 

0.14699 

0.17255 

0.14004 

6.4 

0.09957 

0.13035 

0.14243 

0.14314 

6.6 

0.04352 

0.09700 

0.09680 

0.12888 

6.8 

-0.01592 

0.05073 

0.0411 

0.09857 

7.0 

-0.07104 

-0.00295 

-0.01771 

0.05557 

7.2 

-0.11448 

-0.05745 

-0.07205 

0.00488 

7.4 

-0.14027 

-0.10594 

-0.11471 

-0.04743 

7.6 

-0.14458 

-0.14223 

-0.13988 

-0.09498 

7.8 

-0.12628 

-0.16158 

-0.14384 

-0.13189 

8.0 

-0.08716 

-0.16128 

-0.12556 

-0.15352 
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solutions  are  not  exactly  the  same,  the  order  of  magnitude  and  behavior 
of  results  is  similar.  From  Table  2,  the  same  observations  can  be  made 
for  the  behavior  of  F2.  There  is,  however,  more  error,  quantitatively, 
between  tho  results  for  exact  and  perturbation  methods  and  a  region  of 
qualitative  inaccuracy  between  the  exact  and  perturbation  solutions  exists 
near  t  *  0.  This  takes  the  form  of  a  difference  in  sign  of  Fj  between 
results  from  the  exact  solution  as  compared  to  the  perturbations  solution. 
This  discrepancy  occurred  also  in  the  other  calculations  performed  (not 
shown)  and  will  be  discussed  in  more  detail  later  in  this  chapter. 

In  Tables  3  and  4,  a  comparison  of  results  is  presented  for  modal 
amplitudes  Fj  and  F2  for  a  stable  standing  wave  case.  The  initial 
conditions  for  the  results  in  these  tables  are  Fj^O)  =  0,  F^'tO)  =  1, 

F2(0)  =  0,  F2'(0)  *  0,  G^O)  =  0,  G^O)  *  0,  G2<0)  =  0,  G2'(0)  -  0, 
n  =  40,  c  =  0.1,  and  w  =  0.1.  However,  these  tables  quantitatively 
present  the  effect  of  deviations  of  the  ratio  of  the  second  acoustic 
frequency  to  the  first  from  the  integer  value  of  2  (this  is  controlled 
by  the  parameter  K)  .  These  results  show  that  solutions  for  finite  values 
of  K  are  qualitatively  similar  to  those  for  K  *  0.  This  indicates  that 
the  ratio  of  the  second  acoustic  frequency  to  the  first  does  not  have 
to  be  an  integer  in  order  to  produce  the  type  of  behavior  observed  here. 

A  ratio  near  an  integer  value  will  lead  to  similar  results.  Tables  3 
and  4  also  allow  a  comparison  to  the  results  generated  by  the  program 
in  /ppendix  D  for  the  approximate  analytical  solution  (4.31).  These 
results  presented  in  the  last  column  of  Tables  3  and  4  can  be  compared 
to  the  fourth  column  in  each  of  these  tables  to  determine  the  accuracy 
of  (4.31).  These  comparisons  present  further  evidence  that  the  neglect 
of  gas  dynamic  nonlinearities  does  not  have  an  important  qualitative 
effect. 
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Table  3.  Comparison  of  Results  for  F^  Showing  Effects  of  the  Correction 
Variable  (K)  -  (Ft  =  0,  Fj*  =  1,  F2  *  0,  F2'  =  0,  Gt  *  0,  Qt'  *  0, 

G2  *  0,  G2*  -  0)  -  Stable  Case  (n  *  40)  -  Standing  Waves 


i 

K 

*  1 
=  1 

U  If 

1 

0 

•X  H* 

tl  II 

O  O 

t 

Exact 

Solution 

Perturbation 

Solution 

Exact 

Solution 

Perturbation 

Solution 

Analytic 

Solution 

0.2 

0.19699 

0.18707 

0.19670 

0.19678 

0.19671 

0.4 

0.38335 

0.36396 

0.38172 

0.38210 

0.38186 

0.6 

0.55254 

0.52460 

0.54795 

0.54890 

0.54843 

0.8 

0.69867 

0.66358 

0.68903 

0.69093 

0.69029 

1.0 

0.81667 

0.77635 

0.79957 

0.60302 

0.80234 

1.2 

0.90247 

0.85935 

0.87537 

0.88124 

0.88075 

1.4 

0.95312 

0.91014 

0.91361 

0.92306 

0.92304 

1.6 

0.96699 

0.92744 

0.91310 

0.92740 

0.92820 

1.8 

0.94390 

0.91120 

0.87443 

C. 89468 

0.89666 

2.0 

0.88523 

0.86255 

0.80001 

0.82676 

0.83027 

2.2 

0.79388 

0.78374 

0.69393 

0.72690 

0.73221 

2.4 

0.67411 

0.67805 

0.56163 

0.59954 

0.60682 

2.6 

0.53129 

0.54969 

0.40948 

0.45016 

0.45947 

2.9 

0.37154 

0.40358 

0.24432 

0.28905 

0.29627 

3.0 

0.20133 

0.24520 

0.07300 

0.11102 

0.12386 

3.2 

0.02712 

0.08039 

-0.09782 

-0.06484 

-0.05085 

3.4 

-0.14491 

-0.08439 

-0.26190 

-0.23547 

-0.22097 

3.6 

-0.30901 

-0.24474 

-0.41335 

-0.39410 

-0.37989 

3.8 

-0.45992 

-0.39354 

-0.54662 

-0.53453 

-0.52152 

4.0 

-0.59287 

-0.52617 

-0.65660 

-0.65139 

-0.64055 

4.2 

-0.70357 

-0.63813 

-0.73872 

-0.74028 

-0.73261 

4.4 

-0.78821 

-0.72574 

-0.78928 

-0.79800 

-0.79446 

4.6 

-0.84362 

-0.78623 

-0.80583 

-0.82264 

-0.82407 

4.9 

-0.86746 

-0.81785 

-0.78760 

-0.81364 

-0.82069 

5.0 

-0. 85849 

-0.81988 

-0.73571 

-0.77179 

-0.78490 

5.2 

-0.81682 

-0.79266 

-0.65317 

-0.69920 

-0.71852 

5.4 

-0.74409 

-0.73757 

-0.54464 

-0.59917 

-0.62456 

5.6 

-0.64352 

-0.65697 

-0.41585 

-0.47611 

-0.50704 

5.8 

-0.51967 

-0. 55405 

-0.27311 

-0.33525 

-0.37090 

6.0 

-0.37816 

-0.43279 

-0,12272 

-0.18253 

-0.22172 

6.2 

-0.22516 

-0.29772 

0.02937 

-0.02427 

-0.06554 

6.4 

-0.06694 

-0.15382 

0.17767 

0.13305 

0.09140 

6.6 

0.09056 

-0.0063 

0.31714 

0.28307 

0.24291 

6.8 

0.24192 

0.13957 

0.44298 

0.41977 

0.38308 

7.0 

0.38242 

0.27867 

0.55061 

0,53776 

0.50650 

7.2 

0.50795 

0.40617 

0.63563 

0.63247 

0.60850 

7.4 

0.61493 

0.51774 

0.69408 

0.70029 

0.68527 

7.6 

0.70017 

0.60967 

0 . 72786 

0.73878 

0.73407 

7.8 

0.76090 

0.67898 

0.72020 

0.74667 

0.75327 

8.0 

0.79476 

0.72355 

0.68609 

0.72400 

0.74234 

Table  4.  Comparison  of  Results  for  F2  Showing  Effects  of  the  Correction 
Variable  (K)  -  (Fj  *  0,  Fj'  =  1,  F2  *  0,  F2*  »  0,  Gj  *  0,  Oj'  =  0, 

*  0,  62'  *  0)  -  Stable  Case  (n  *  40)  -  Standing  Waves 


i  * 
K  * 

1 

1 

i  a  3 
K  =  C 

i  *  0 

K  s  0 

t 

Exact 

Solution 

Perturbation 

Solution 

Exact 

Solution 

Perturbation 

Solution 

Analytic 

Solution 

0.2 

0.00011 

-0.00400 

0.00015 

-0.00419 

-0.00192 

0.4 

0.00099 

-0.00996 

0.00136 

-0.01035 

-0.00696 

0.6 

0.00354 

-0.01599 

0.00479 

-0.01597 

-0.01337 

0.8 

0.00860 

-0.02013 

0.01141 

-0.01856 

-0.01885 

1.0 

0.01665 

-0.02071 

0.02149 

-0.01626 

-0.02113 

1.2 

0.02761 

-0.01669 

0.03443 

-0.00838 

-0.01856 

1.4 

0.04072 

-0.00793 

0.04866 

0.00438 

-0.01059 

1.6 

0.05458 

0.00479 

0.06180 

0.01995 

0.00208 

1.8 

0.06726 

0.01983 

0.07113 

0.03529 

0.01749 

2.0 

0.07663 

0.03491 

0.07416 

0.04696 

0.03278 

2.2 

0.08072 

0.04750 

0.06910 

0.05185 

0.04472 

2.4 

0.07803 

0.05521 

0.05549 

0.04793 

0.05038 

2.6 

0.06793 

0.05622 

0.03439 

0.03468 

0.04776 

2.8 

0.05081 

0.04957 

0.00841 

0.01343 

0.03626 

3.0 

0.02816 

0.03542 

-0.01367 

-0.01282 

0.01697 

3.2 

0.00244 

0.01507 

-0.04246 

-0.03978 

-0.00745 

3.4 

-0.02323 

-0.00917 

-0.05871 

-0.06265 

-0.03313 

3.6 

-0.04544 

-0.03425 

-0.06410 

-0.07695 

-0.05561 

3.8 

-0.06099 

-0.05679 

-0.05687 

-0.07938 

-0.07067 

4.0 

-0.06735 

-0.07356 

-0.03727 

-0.06854 

-0.07507 

4.2 

-0.06307 

-0.08196 

-0.00764 

-0.04530 

-0.06722 

4.4 

-0.04302 

-0.08040 

0.02778 

-0.01278 

-0.04759 

4.6 

-0.02355 

-0.06854 

0.06351 

0.02408 

-0.01872 

4.8 

0.00767 

-0.04745 

0.09363 

0.05926 

0.01509 

5.0 

0.04187 

-0.01945 

0.11279 

0.08668 

0.04845 

5.2 

0.07470 

0.01210 

0.11719 

0.10127 

0.07573 

5.4 

0.10173 

0.04330 

0.10529 

0.09987 

0.09206 

5.6 

0.11915 

0.07013 

0.07824 

0.08191 

0.09414 

5.8 

0.12425 

0.08905 

0.03976 

0.04959 

0.08095 

6.0 

0.11586 

0.09745 

-0.00434 

0.00764 

0.05395 

6.2 

0.09462 

0.09402 

-0.04711 

-0.03736 

0.01700 

6.4 

0.06286 

0.07894 

-0.08154 

-0.07807 

-0.02424 

6.6 

0.02438 

0.05389 

-0.10178 

-0.10757 

-0.06318 

6.8 

-0.01630 

0.02182 

-0.10409 

-0.12058 

-0.09336 

7.0 

-0.05322 

-0.01335 

-0.08754 

-0.11433 

-0.10952 

7.2 

-0.08237 

-0.04728 

-0.05423 

-0.08921 

-0.10854 

7.4 

-0.09956 

-0.07576 

-0.00908 

-0.04879 

-0.09001 

7.6 

-0.10237 

-0.09523 

0.04092 

0.00093 

-0.05639 

7.8 

-0.09018 

-0.10322 

0.08777 

0.05190 

-0.01267 

8.0 

-0.06428 

-0.09872 

0.12375 

0.09586 

0.03434 

In  Tables  5  and  6,  a  comparison  of  results  are  presented  for  modal 
amplitudes  Fj  and  Fj  for  an  unstable  standing  wave  showing  the  effect  of 
neglecting  gas-dynamic  nonlinearities.  It  can  be  seen  that  the  gas 
dynamic  nonlinearities  have  little  qualitative  effect  on  the  results. 

In  Tables  7  and  8,  a  comparison  of  results  are  presented  for  modal 
amplitudes  Fj  and  Fj  for  an  unstable  standing  wave  case  showing  the  effects 
of  K.  The  results  for  zero  and  non-zero  are  qualitatively  similar. 

These  tables  are  representative  of  the  cases  that  were  investigated 
in  the  course  of  this  research.  Only  cases  involving  standing  waves  were 
presented.  The  same  behavior,  however,  can  be  observed  for  the  cases 
involving  traveling  waves. 

In  Table  9,  a  comparison  of  stability  boundaries  is  presented 
based  upon  the  interaction  index  (n)  which  is  a  measure  of  the  strength 
of  the  combustion  process.  For  standing  waves  and  the  given  conditions 
shown,  the  stability  limit  for  a  process  with  gas  dynamic  nonlinearities 
considered  and  K  =  0  is  between  45-50.  When  both  gas  dynamic  non- 
linearities  and  the  correction  variable  are  considered,  the  stability 
limit  is  increased  to  67.5-69.  Finally,  when  considering  only  the 
correction  variable  with  no  gas-dynamic  non-linearity  effect,  the  stability 
limit  is  72-72.5.  The  results  show  that  the  neglect  of  gas  dynamic 
non  linearities  slightly  underestimates  the  stability  boundary  and  that 
the  increasing  K  increases  the  stability  limit. 

In  Table  10,  a  comparison  of  stability  boundaries  is  presented 
based  upon  the  interaction  index  for  traveling  waves.  These  results  provide 
additional  confirmation  of  the  conclusions  discussed  in  the  previous 
paragraph  and  also  illustrate  the  fact  that  standing  waves  are  roughly 
twice  as  stable  as  traveling  waves.  This  is  consistent  with  the 
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Table  5.  Comparison  of  Results  for  Fj  Showing  Effect  of  the  Gas  Dynamic 
Index  (i)  -  (Fj  »  0,  Fj'  *  1,  Fj  8  0,  Fj'  8  0,  Gj  8  0.  G^'  8  0,  G2  8  0, 
G2'  8  0)  -  Unstable  Case  (n  =  75)  -  Standing  Waves 


H  H 

II  II 

•H 

i  8  0 

K  =  1 

t 

Exact 

Perturbation 

Exact 

Perturbation 

Solution 

Solution 

Solution 

Solution 

0.2 

0.19699 

0.18717 

0.19699 

0.18703 

0.4 

0.38335 

0.38215 

0.38336 

0.36402 

0.6 

0.55250 

0.52613 

0.55258 

0.52513 

0.8 

0.69847 

0.69244 

0.69881 

0.66518 

1.0 

0.81592 

0.80625 

0.81700 

0.77974 

1.2 

0.90030 

0.86736 

0.90294 

0.86535 

1.4 

0.94803 

0.92126 

0.95337 

0.91952 

1.6 

0.95671 

0.94184 

0.96695 

0.94089 

1.8 

0.92557 

0.92873 

0.94003 

0.92921 

2.0 

0.85572 

0.88276 

0.87584 

0.88529 

2.2 

0.75039 

0.80579 

0.77580 

0.81106 

2.4 

0.61484 

0.70079 

0.64411 

0.70941 

2.6 

0.45593 

0.57161 

0.48677 

0.58410 

2.8 

0.28147 

0.42287 

0.31118 

0.43963 

3.0 

0.09945 

0.25987 

0.12554 

0.28107 

3.2 

-0.08272 

0.08829 

-0.06194 

0.11388 

3.4 

-0.25864 

-0.08588 

-0.24371 

-0.05620 

3.6 

-0.42307 

-0.25669 

-0.41342 

-0.22350 

3.8 

-0.57177 

-0.41827 

-0.56608 

-0.38241 

4.0 

-0.70107 

-0.56512 

-0.69790 

-0.52773 

4.2 

-0.80745 

-0.69223 

-0.80594 

-0.63473 

4.4 

-0.88721 

-0.79547 

-0.88767 

-0.75934 

4.6 

-0.93641 

-0.87123 

-0.94056 

-0.83825 

4.8 

-0.95124 

-0.91703 

-0.96203 

-0.88902 

5.0 

-0.92862 

-0.93138 

-0.94961 

-0.91020 

5.2 

-0.86710 

-0.91381 

-0.90146 

-0.90121 

5.4 

-0.76771 

-0.86492 

-0.81711 

-0.86255 

5.6 

-0.63434 

-0.78635 

-0.69815 

-0.79563 

5.8 

-0.47368 

-0.68071 

-0.54872 

-0.70274 

6.0 

-0.29439 

-0.55149 

-0.37551 

-0.58703 

6.2 

-0.10598 

-0.40299 

-0.18717 

-0.45231 

6.4 

0.08258 

-0.22558 

0.00672 

-0.30299 

6.6 

0.26374 

-0.06015 

0.19691 

-0.14391 

6.8 

0.43198 

0.10715 

0.37555 

0.01979 

7.0 

0.58355 

0.28009 

0.53682 

0.18290 

7.2 

0.71613 

0.44478 

0.67707 

0.34025 

7.4 

0.82779 

0.59627 

0.79428 

0.48686 

7.6 

0.91635 

0.72909 

0.88737 

0.61813 

7.8 

0.97878 

0.83897 

0.95535 

0.72996 

8.0 

1.01111 

0.92224 

0.99664 

0.81891 

Table  6.  Comparison  of  Results  for  F2  Showing  Effect  of  the  Gas 
Dynamic  Index  (i)  -  (Fj  =  0,  Fj'  s  1,  Fj  s  0,  F2'  =  0,  Gj  =  0, 

Gi'  a  0,  G2  =  0,  G2'  =  0)  -  Unstable  Case  -  (n  =  75)  -  Standing  Waves 


H  H 

ft  II 

O 

II  II 

t 

Exact 

Perturbation 

Exact 

Perturbation 

Solution 

Solution 

Solution 

Solution 

0.2 

0.00013 

-0.00548 

0.00003 

-0.00316 

0.4 

0.00124 

-0.01544 

0.00054 

-0.01172 

0.6 

0.00473 

-0.02691 

0.00257 

-0.02333 

0.8 

0.01210 

-0.03643 

0.00752 

-0.03484 

1.0 

0.02443 

-0.04077 

0.01670 

-0.04288  j 

1.2 

0.04197 

-0.03749 

0.03089 

-0.04728 

1.4 

0.06387 

-0.025518 

0.04994 

-0.03793 

1.6 

0.08806 

-0.005425 

0.07263 

-0.02258 

1.8 

0.11143 

0.02061 

0.0966 

0.00033 

2.0 

0.13027 

0.049007 

0.1186 

0.02813 

2.2 

0.14082 

0.07537 

0.13504 

0.05696 

2.4 

0.14005 

0.09513 

0.14236 

0.00237 

2.6 

0.12628 

0.10426 

0.13796 

0.09995 

2.8 

0.09967 

0.10004 

0.12064 

0.10606 

3.0 

0.06242 

0.08152 

0.09109 

0.09851 

3.2 

0.01860 

0.04984 

0.05192 

0.07688 

3.4 

-0.02636 

0.008166 

0.0076 

0.04283 

3.6 

-0.06636 

-0.03858 

-0.03628 

-0.00096 

3.8 

-0.09548 

-0.08436 

-0.07362 

-0.046816 

4.0 

-0.10889 

-0.12285 

-0.09877 

-0.09134 

4.2 

-0.10354 

-0.14826 

-0.10733 

-0.12758 

4.4 

-0.07873 

-0.15626 

-0.09688 

-0.15016 

4.6 

-0.03636 

-0.14461 

-0.06742 

-0.15524 

4.8 

0.01909 

-0.11361 

-0.02153 

-0.14104 

5.0 

Q. 08104 

-0.06617 

0.03583 

-0.10827 

5.2 

0.14148 

-0.00758 

0.09778 

-0.06010 

5.4 

0.19206 

0.055118 

0.15638 

-0.00186  | 

5.6 

0.22528 

0.11399 

0.20358 

0.05954 

5.8 

0.23565 

0.16121 

0.23235 

0.11651 

6.0 

0.22067 

0.19016 

0.23779 

0.16168 

6.2 

0.18130 

0.19633 

0.21796 

0.18888 

6.4 

0.12201 

0.17793 

0.17434 

0.19397 

6.6 

0.05019 

0.13633 

0.11182 

0.17547 

6.8 

-0.02482 

0.07595 

0.03808 

0.13477 

7.0 

-0.09302 

-0.02526 

-0.03737 

0.07614 

7.2 

-0.14503 

-0.07144 

-0.10448 

0.006209 

7.4 

-0.17338 

-0.140235 

-0.15407 

-0.06680 

7.6 

-0.17347 

-0.19368 

-0.17900 

-0.13401 

7.8 

-0.14428 

-0.22450 

-0.17519 

-0.18700 

8.0 

-0.08859 

-0.22808 

-0.14216 

-0.21888 
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Table  7.  Comparison  of  Results  for  Fj  Showing  the  Effects  of  the 
Correction  Variable  (K)  -  (Fj  =  0,  F^'  =  1,  F2  =  0,  F2*  8  0,  Gj  =  0, 
Gj/  *  0,  G2  8  0,  G2'  8  0)  -  Unstable  Cases  (n  *  70)  -  Standing  Waves 


II  II 

1 

1 

i  8 
K  = 

1 

0 

i  =  0 

K  8  0 

t 

Exact 

Perturbation 

Exact 

Perturbation 

Analytic 

Solution 

Solution 

Solution 

Solution 

Solution 

0.2 

0.19699 

0.187155 

0.19670 

0.19687 

0.19675 

0.4 

0.38335 

0.36462 

0.38172 

0.38261 

0.38217 

0.6 

0.55251 

0.52587 

0.54791 

0.55028 

0.54942 

0.8 

0.69850 

0.66615 

0.68878 

0.69371 

0.69249 

1.0 

0.81604 

0.78072 

0.79865 

0.80767 

0.80630 

1.2 

0.90066 

0.86596 

0.87280 

0.88811 

0.88696 

1.4 

0.94887 

0.91928 

0.90775 

0.93226 

0.93184 

0.95842 

0.93925 

0.90166 

0.93878 

0.93468 

0.92863 

0.92554 

0.85477 

0.90772 

0.91061 

0.86067 

0.87904 

0.76962 

0.84061 

0.84612 

0.75771 

0.80166 

0.65109 

0.74802 

0.74903 

0.62481 

0.69648 

0.50590 

0.61107 

0.62333 

0.46860 

0.56735 

0.34194 

0.45807 

0.47407 

2.8 

0.29658 

0.41903 

0.16733 

0.26849 

0.30714 

3.0 

0.11649 

0.25679 

-0.01027 

0.10616 

0.12906 

3.2 

-0.06441 

0.08639 

-0.18416 

-0.05313 

-0.05327 

3.4 

-0.23973 

-0.08618 

-0.34860 

-0.25962 

-0.23280 

3.6 

-0.40414 

-0.25497 

-0.49854 

-0.42949 

-0.40262 

3.8 

-0.55318 

-0.41415 

-0.62912 

-0.58151 

-0.55621 

4.0 

-0.68301 

-0.55829 

-0.73533 

-0.70960 

-0.68766 

4.2 

-0.78998 

-0.68250 

-0.81197 

-0.80859 

-0.79192 

4.4 

-0.87041 

-0.78258 

-0.85399 

-0.87449 

-0.86494 

4.6 

-0.92055 

-0.85525 

-0.85736 

-0.90458 

-0.90388 

4.8 

-0.93689 

-0.89811 

-0.82008 

-0.89755 

-0.90718 

-0.91671 

-0.90981 

-0.74306 

-0.85356 

-0.87463 

-0.85886 

-0.89005 

-0.63045 

-0.80213 

-0.80736 

-0.76447 

-0.83962 

-0.48925 

-0.66264 

-0.70785 

-0.63724 

-0.76029 

-0.32818 

-0.52309 

-0.57982 

5.8  I 

-0.48340 

-0.65486 

-0.15634 

-0.36111 

-0.42807 

-0.31100 

-0.52691 

0.01809 

-0.18308 

-0.25835 

-0.12888 

-0.38083 

0.18852 

-0.07645 

-0.07712 

0.05445 

-0.22156 

0.34995 

0.10764 

0.10864 

6.6 

0.23170 

-0.05445 

0.49844 

0.37422 

0.29176 

6.8 

0.39723 

0.11485 

0.63023 

0.46032 

0.46509 

7.0 

0.54707 

0.28072 

0.74089 

0.69003 

0.62177 

7.2 

0.67846 

0.43759 

0.82500 

0.81052 

0.75551 

7.4 

0.78917 

0.58027 

0.87632 

0.89873 

0.86085 

7.6 

0.87687 

0.70394 

0.88888 

0.95061 

0.93330 

7.8 

0.93868 

0.80457 

0.85847 

0.96357 

0.96962 

8.0 

0.97109 

0.87877 

0.78419 

0.93647 

0.96786 

I  S-  C0"P"'1"0"  of  R“ul«  to  F2  Showing  th.  Effects  of  the 
Correction  Vsriable  (K)  -  (F,  .  0,  r/  .  1.  p,  .  o,  r2> 

Si  0,  Cj  .  0,  8j  .  0)  -  Unstable  Cases  <„  .  70)  .  Sta„dlng  ^ 


0.2 

0.1* 

0.6 

0.8 

1.0 

1.2 

1.4 
1.6 
1.8 
2.0 
2.2 

2.4 
2.6 
2.8 
3.0 

3.2 

3.4 

3.6 

3.8 
4.0 

4.2 

4.4 

4.6 

4.8 
5.0 

5.2 

5.4 

5.6 

5.8 
6.0 

6.2 

6.4 

6.6 

6.8 
7.0 
7.2 

7.4 
7.6 
7.8 
8.0 


Exact 

Solution 

0.00013 
0.00120 
0.00456 
0.01160 
0.02332 
0.03992 
0.06057 
0.08330 
!  0.10517 

0.12271 
0.13242 
0.13150 
0.11839 
0.09328 
0.05822 
0.01701 
-0.02528 
-0.06295 
-0.09049 
-0.10337 
-0.09872 
-0.07590 
-0.03665 
0.01487 
0.07254 
0.12595 
0.17637 
0.20783 
0.21824 
0,20519 
0.16949 
0.11514 
0.04883 
-0.02090 
-0.08483 
-0.13429 
-0.16228 
-0.16442 
-0.13960 
-0.09020 


Perturbation 

Solution 


Exact 

Solution 


Perturbation 

Solution 


i  =  0 
K  «  0 

Analytic 

Solution 


-0.00527 

-0.01466 

-0.02534 

-0.03411 

-0.03789 

-0,03451 

-0.02299 

-0.00396 

0.02047 

0.04692 

0.07126 

0.08922 

0.097153 

0.11651 

0.07467 

0.04468 

0.00567 

-0.03775 

-0.07994 

-0.11504 

-0.13779 

-0.1442 

-0.13261 

-0.1032 

-0.05894 

-0.00478 

0.05266 

0.10617 

0.14862 

0.17409 

9.17865 

0.16089 

0.12229 

0.06705 

0.00165 

-0.06594 

-0.12723 

-0.17433 

-0.20087 

-0.21857 


0.00017 
0.00165 
0.00619 
0.01544 
0.03024 
0.05013 
0.07306 
0.09555 
0.11324 
0.12170 
0.11751 
0.09917 
0.06770 
0.02679 
-0.01761 
-0.05827 
-0.08789  | 
-0.10041  I 
-0.09221 
-0.06290 
-0.01566 
0.04299 
0.10404 
0.15731 
0.19315 
0.20424 
0.18706 
0.14264 
0.07764 
1  0.00147 
-0.0733 
-0.13397 
-0.16972 
-0.17346 
-0.14308 
-C. 08205 
0.00095 
0.09312 
0.17949 
0.24523 


-0.00562 

-0.01557 

-0.02603 

-0.03276 

-0.03218 

-0.02232 

-0.01862 

0.021911 

0.04915 

0.07264 

0.08675 

•0.087046 

0.071414 

0.04063 

-0.001408 

-0.04828 

-0.09196 

-0.12426 

-0.13827 

-0.12986 

-0.09854 

-0.04785 

0.01493 

0.07988 

0.13602 

0.17311 

0.1Q3575 

0.163922 

0.11562 

0.04509 

-0.03703 

-0.11757 

-0.182899 

-0.22122 

-0.22468 

-0.19089 

-0.12367 

-0.03255 

0.06839 

0.16276 


-0.00336 
-0.01219 
-0.02343 
-0.03305 
-0.03709 
-0.03264 
-0.01865 
0.00367 
0.03093 
0.05810 
0.07947 
0.08978 
0.08535 
0.06502 
0.03053 
-0.01346 
-0.06005 
-0.10122 
-0.12919 
-0.13785 
-0.12403 
-0.08825 
-C. 03490 
-0.02828 
0.09131 
0.14357 
0.17558 
0.18070 
0.15640 
0.10495  j 
0.03331  I 
-0.04784 
!  -0.12564 
-0.18709 
-0.22125 
-0.22111 
-0.18495 
-0.11690 
-0.02652 
0.07753 
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T«bl«  9.  Comparison  of  Stability  Boundarias  Based  on  tha  Interaction 
Index  (n)  -  (F*  *  0,  Fj*  ■  1,  F2  *  0,  Fj*  *  0,  Qj  *  0»  G*’  ■  0, 

C2  *  0,  G2’  *  0)  -  Standing  Waves  -  Epsilon  -  0.1 


Gas  Dynamic  Index 
Correction  Variable 

Stability  Boundaries  "j 

Exact  Solution 
n  -  Stable  -  Unstable 

Perturbation  Solution 
n  -  Stable  -  Unstable 

i  «  1 
K  «  1 


67.5  -  69 


67.5  -  69 
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Table  10.  Comparison  of  Stability  Boundaries  Based  on  the  Interaction 
Index  (n)  -  (Fj  »  0,  Fj*  *  -1,  F2  »  0,  F2'  *  0,  Gj  *  1,  Gj*  «  0, 

G2  *  0,  G2'  a  0)  -  Traveling  Waves  -  Epsilon  -  0.1 


Gas  Dynamic  Index 
Correction  Variable 

Stability  Boundaries 

Exact  Solution 
n  -  Stable  -  Unstable 

Perturbation  Solution 
n  -  Stable  -  Unstable 

i  a  1 

K  a  1 

27.5  -  28 

31.5  -  32 

i  =  0 

K  a  1 

30  -  31 

36.35  -  36.5 

i  a  1 

K  =  0 

25  -  30 

25-30 
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approximate  analytical  stability  aquations  (4.31)  and  (4.52).  The 
perturbation  method  tends  to  predict  slightly  higher  stability  limits 
than  the  exact  solution  method  for  both  standing  and  traveling  waves. 
Within  the  accuracy  of  the  tabulated  values,  this  is  apparent  only  in 
the  first  two  rows  of  Table  10. 

In  Table  11,  a  comparison  of  the  effect  of  different  initial 
conditions  imposed  on  the  stability  boundaries  for  both  standing  and 
traveling  waves  is  presented.  From  the  results  of  two  sets  of  initial 
conditions  for  each  case,  it  can  be  seen  that  the  varying  of  initial 
conditions  has  no  significant  effect  on  the  stability  boundaries  for 
both  standing  waves  or  traveling  waves. 

In  Table  12,  the  variation  of  the  stability  limit  with  e  is 
presented  for  standing  waves.  From  Table  12,  the  results  show  that  the 
smaller  the  term  epsilon  the  greater  the  stability  limit.  Therefore, 
the  order  term  has  a  significant  effect  on  the  interaction  index.  In 
Chapter  4,  a  relation  was  proposed  for  the  case  of  i  *  0  and  K  =  0 
which  was  n  =  C/e  where  C  is  a  constant.  Assuming  the  validity  of  the 
relation,  the  values  for  this  constant  are  given  for  each  given  epsilon 
and  interaction  index.  This  shows  that,  in  general,  C  is  a  weak  function 
of  e. 

In  Table  13,  a  comparison  of  the  effect  of  e  is  presented  for 
traveling  waves  when  both  gas  dynamic  nonlinearities  and  correction 
variables  are  considered.  Again,  the  results  show  that  the  smaller  the 
term  epsilon,  the  greater  the  stability  limit.  The  perturbation  method 
again  predicts  slightly  greater  stability  limits  than  does  the  exact 
solution  method.  Therefore,  again,  the  order  term  has  a  strong  effect 
concerning  the  stability  of  combustion. 
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Table  11.  Comparison  of  ths  Effect  of  Different  Initial  Conditions 
Imposed  for  Standing  and  Traveling  Waves  for  i  «  1  and  K  *  1 
Epsilon  *  0.1 

(a)  Standing  Waves  -  1.  Fj  *  0,  '  *  1»  F2  *  0,  F2’  *  0 

Gj  *  0,  •  0,  G2  «  0,  G2*  »  0 


2.  F*  a  1,  Fj'  *  0,  F2  =  0,  F2’  «  0 
Gj  8  0,  Gj*  g  0,  G2  «  0,  G2 1  «  0 


Initial  Condition 

Sets 

Stability  Boundaries 

Exact  Solution 
n  -  Stable  -  Unstable 

Perturbation  Solution 
n  -  Stable  -  Unstable 

1. 

67.5  -  69 

67.5  -  69 

2. 

65  -  70 

65  -  70 

(b)  Traveling  Waves  -  1.  F^  =  0,  Fj/  =  -1,  F2  =  0,  F2'  =  0 

G1  *  1*  Gi*  s  0,  G2  =  0,  G2'  =  0 


2.  F*  =  1,  F1’  a  0,  F2  a  0,  F2*  a  0 
Gj  »  0,  *  *  "1 1  G2  *  0 ,  G2  *  0 


Initial  Condition 

Sets 

Stability  Boundaries 

Exact  Solution 
n  -  Stable  -  Unstable 

Perturbation  Solution 
n  -  Stable  -  Unstable 

1. 

27.5  -  28 

31.5  -  32 

2. 

27.5  -  28.5 

31  -  31.5 

Table  13.  Comparison  of  the  Effects  of  the  Order  Term  Epsilon  (F 
F2  ~  °»  G!  -  1*  Gj  =  0,  G2  =  o,  G  ’  =  0)  Traveling  Waves  -  when 
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Thus,  from  these  representative  tables  of  results ,  it  is  observed 
that  the  correction  variable  is  important  in  the  stability  of  standing 
waves,  but  does  not  play  a  major  role  in  the  stability  of  traveling  waves. 
It  is  observed  that  the  gas  dynamic  nonlinearities  seem  to  have  little 
influence  on  the  stability  of  either  standing  or  traveling  waves.  It 
is  observed  that  initial  conditions  of  the  modal  amplitudes  have  little 
or  no  influence  in  the  stability  of  either  standing  or  traveling  waves. 

And  finally,  it  is  observed  that  the  order  term  epsilon  and.  the  inter¬ 
action  index  governing  the  strength  of  combustion  in  the  process  are 
strongly  coupled  thus  affecting  the  limits  of  stability. 

Before  completing  this  chapter,  it  is  desired  to  investigate  the 
sign  discrepency  mentioned  previously  between  the  exact  and  perturbation 
solutions  for  f  which  occur  near  t  =  0.  For  simplicity,  it  will  be 
assumed  that  i  *  K  *  0  and  that  for  t  <<  1  the  first  modal  amplitude  can 
be  represented  with  sufficient  accuracy  by  fj  =  sint.  Then,  the 
equation  for  f2  will  be  solved  and  the  result  simplified  for  t  <<  1. 

This  will  be  done  first  for  w  =  0  and  then  for  w  t  0.  For  «  «  0, 

(3.21)  leads  to 


d2f 

dt^  *  4f2 


»  k  cun 


cos2t 


(5.1) 


with  initial  conditions 


f2(0)  =  0 


f2  (0)  *  0. 
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Evalu  i.h*  homogeneous  and  particular  solutions  by  tha  usual  mannar 
and  avaluating  tha  constants,  tha  results  become 


1__ 

16 


cun 


t- 


cos2t  -  t  sin2t 


l 


(5.2) 


In  terms  of  the  perturbation  parameters  (4.1),  equation  (5.2)  can  be 
written  as 


1_ 

16 


un  e(l  -  cos2C)  -  n  sin2£ 


3. 


(5.3) 


To  the  order  of  approximation  e  which  the  perturbation  solution  should 
model,  equation  (5.3)  becomes 

f2  =  *  16  “nn  Sin25  +  0(e)‘  (5.4) 

By  expanding  equation  (5.2)  into  a  Taylor  series  expansion  of  three  terms, 
equation  (5.2)  becomes 


f 


2 


24 


-  4 

eunt  + 


(5.5) 


which  is  always  positive. 

Therefore,  tha  exact  method  for  small  time  will  yield  fj  modal 
amplitude  always  as  a  positive  quantity. 

By  imposing  identical  conditions  to  the  perturbation  equations 


(4.12),  the  result  becomes 
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uin 


(5.6) 


with  the  condition 


b2(o)  *  0. 

Solving  equation  (5.6), 

B2  =  ’  IF  “nn*  (5.7) 


Recalling  that  f2  a  B2  sin2£,  the  result  becomes 

f2  *  “  IF  5nn  sin2?  +  0(c)  (5.8 

which  is  identical  to  the  result  of  equation  (5.4)  for  the  wave  equation 
solution.  Thus,  the  perturbation  method  gives  the  correct  result.  It 
can  be  seen  that  for  t  «  1  the  exact  solution  predicts  a  positive  f2 
and  by  inspection  of  equation  (5.8),  the  perturbation  method  predicts  a 
negative  f^.  This  is  precisely  the  behavior  observed  in  the  numerical 
solutions. 

For  u  /  0,  a  similar  analysis  can  be  performed.  The  appropriate 
equation  for  f2  is  now 

d2f  df  r 

+  u  +  4f2  =  ,CE‘*’nj  1  "  cos2t 


(5.9) 


”,  "rvTrrWFy 
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with  conditions 


f2(0)  *  0 
f2'(0)  *  0. 


Solving  the  homogeneous  and  particular  solution  by  the  usual  manner  and 
evaluating  the  appropriate  constants  the  result  becomes 


f2  *  e 


-w/2t 


[-  rs eSn  cos 


en  re  -  a2  i 
16  L/i6  -_a2  J 


/l6  -  5^  1  - 

sxn  —————  t  +  ?r  ewn 
«  16 


en  sin2t 


’] 


(5.10) 


Expanding  (5.10)  for  small  co  into  the  appropriate  Taylor  series,  expanding 
and  neglecting  terms  of  0(w)  leads  to 


f2  * 


fr  cSn  [; 


1  -  cos2t  -  t  sin2tj 


(5.11) 


which  is  identical  to  (5.2). 

By  imposing  the  identical  conditions  on  the  perturbation  equation 
(4.12),  the  resulting  equation  become 

dB 

3^-  +  h  oB2  *  -  jg-  un  (5.12) 


with  the  condition 


B2(0)  =  0. 
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Solving  equation  (5.12)  by  the  usual  manner,  evaluating  the  constants, 
and  transforming  the  perturbation  variables  to  real  time  variables 

f2  =  '  —  [l  -  e  sin  2t  .  (5.13) 

This  is  always  negative  for  t  «  1.  Expanding  the  exponential  function 
by  the  Taylor  series  expansion  and  neglect  terms  of  o(ui)  leads  to 

f2  =  n  sin  2£  +  0(e)  (5.14) 

which  is  identical  to  (5.8). 

To  observe  the  behavior  of  equation  (5.10)  for  small  time, 
expand  this  equation  into  a  Taylor  series  of  0 ( t14 ) .  Expanding  and 
grouping  terms  according  to  their  order  of  magnitude,  the  terms  of 
0(1),  0(t),  0(t2),  0(t3)  vanish.  Therefore,  f2  is  comprised  of  terms 
from  Oft1*)  which  is 


f  =  en£St4  L  +  3^1  +  “Ill 

r2  24  L  8  64  J  .  (5.15) 

4 

Again,  for  any  small  time  t,  fj  is  always  positive  since  t  is  always 
positive.  Neglecting  higher  powers  of  w,  the  resulting  equation  becomes 
equation  (5.5)  for  the  undamped  case.  Again  it  can  be  seen  that  the 
exact  and  perturbation  methods  predict  opposite  signs  for  f2  when  t  «  1. 
These  results  are  based  on  approximations  and  cannot  be  considered 


definitive.  They  do,  however,  lend  plausibility  to  the  numerical  results 
discussed  earlier.  It  is  believed  that  this  sign  discrepancy  is  due 
to  the  inability  of  the  perturbation  solution  to  accurately  represent 
the  exact  solution  for  t  «  1  and  not  due  to  any  error  in  the  computer 
program  used  to  compute  the  perturbation  solution. 
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Chapter  6 

CONCLUSION  AND  RECOMMENDATIONS 

The  primary  objective  of  this  presentation  has  been  the  development 
of  analytical  techniques  to  solve  the  problem  of  combustion  instabilities 
occurring  in  an  annular  combustion  chamber.  The  analytical  techniques 
used  were  the  modified  Galerkin  method  applied  to  the  acoustic  wave 
equations  which  yielded  a  set  of  time-dependent  modal  amplitude  equations 
and  the  two-variable  perturbation  method  which  yield  a  set  of  time- 
dependent  equations  which  approximated  the  behavior  of  the  first  set  of 
equations.  Both  methods  produced  results  which  were  relatively  easy  to 
apply  and  used  the  Runge-Kutta  algorithm  which  required  little  computation 
time.  An  alternative  approach  to  solve  this  problem  would  be  a  finite 
difference  approach.  However,  difficulties  can  be  foreseen  in  the 
development  of  the  finite  difference  equations  modelling  the  problem 
along  with  the  complications  occurring  due  to  the  boundary  conditions  of 
the  problem.  Thus,  the  benefits  of  the  methods  discussed  in  this  thesis 
can  be  appreciated. 

From  the  numerical  and  graphical  presentation  of  results  in  Chapter 
5,  the  following  observations  can  be  made.  First,  the  effect  of  the  gas- 
dynamic  nonlinearities  seems  to  be  small  in  both  methods  of  analysis  for 
velocity  sensitive  combustion.  This  point  can  be  observed  from  a 
quantitative  comparison  of  the  tabular  results  or  by  observing  the  effects 
of  this  condition  on  the  stability  boundaries.  Second,  the  effect  of  the 
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correction  variable  modelling  the  physical  boundaries  of  the  chamber  seems 
to  have  a  significant  effect  in  both  methods  of  analysis  for  velocity 
sensitive  combustion.  By  including  the  effect  of  this  correction  variable, 
a  significant  increase  occurs  in  the  interaction  index  which  is  the 
criteria  for  the  stability  of  the  system.  However,  this  effect  seems  to 
be  more  significant  for  the  standing  wave  case  than  the  traveling  wave 
cases.  The  effects  of  initial  conditions  for  the  time  dependent  equations, 
the  numerical  value  for  the  burning  rate  and  step  size  of  integration , 
seem  to  have  very  little  significance  in  the  measure  of  the  stability 
limits  of  velocity  sensitive  combustion.  However,  the  order  term  epsilon 
has  a  strong  effect  upon  the  stability  of  the  problem.  This  is  to  be 
excepted  since  the  order  term  is  the  measure  of  the  effect  of  non- 
linearities  occurring  in  the  system.  The  increase  in  this  value  corresponds 
to  a  decrease  in  the  stability  limit  which  is  physically  reasonable. 

In  this  study,  the  effect  of  time  delay  of  the  combustion  process 
was  neglected.  However,  time  delay  has  been  found  in  other  studies  to 
be  an  important  phenomena  in  correctly  modelling  the  actual  problems  of 
velocity  sensitive  combustion.  It  is  recommended  that  this  effect  can 
be  incorporated  by  including  the  corresponding  terms  with  j  *  1  in  the 
acoustic  wave  equations  (3.20).  A  corresponding  set  of  perturbations  can 
then  be  derived  to  account  for  time  delay  and  both  these  equations  and 
equations  (3.20)  can  be  numerically  evaluated  by  modifing  the  existing 
Runge-Kutta  programs  presented  in  the  Appendices.  It  is  also  recommended 
that  an  experimental  program  be  developed  to  measure  the  effects  of 
velocity  sensitive  combustion  in  an  annular  combustion  chamber.  Once 
achieving  this  goal,  one  could  correlate  the  measurement  results  to  the 
analytical  results  that  have  been  presented  to  ascertain  the  validity  of 
this  analysis. 
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Since  instability  of  combustion  is  sensitive  to  small  changes  in 
engine  geometry  and  operating  conditions,  a  particular  engine  must  be 
subjected  to  a  large  number  of  firings  before  its  designers  can  say 
confidently  that  it  is  free  from  instability.  With  a  large  engine  such 
testing  can  account  for  a  substantial  part  of  development  costs.  Herein 
lies  the  importance  of  devising  reliable  theories  of  instability  and 
inexpensive  tests  of  a  propellant's  acoustical  characteristics.  Until 
instability  of  combustion  is  understood  well  enough  so  that  it  can  be 
eliminated  while  an  engine  is  in  the  design  stage,  rocket  engines  must 
continue  to  be  intensively  tested  for  stability— particularly  when 
the  lives  of  astronauts  will  eventually  depend  on  safe,  reliable 
operation  of  the  engine  [17]. 
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GENERAL  TIME  DELAY  FUNCTION 

The  development  end  nature  of  the  time-delay  function  is  of  the 
same  form  of  the  convolution  integral  for  impulse  response  in  vibration 
theory.  The  general  form  of  the  time  delay  function  is 


/  d“o 

w(t)  ■  I  J(t  -  O  dc  . 


(A.l) 


A  simple  illustration  of  the  time  delay  function  is  in  the  case  of  a 

finite  step  function  J(t). 

J(t)| 


(some  specific  time  constant) 
Figure  Al.  Step  Function  J(t) 

From  the  figure,  the  step  function  J(t)  is  defined  as 


J(t)  = 


1  t  <  T 


0  t  >  T 


(A. 2) 


Therefore,  substituting  some  time  delay  (t  -  C)  for  time  t,  the  result  is 
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Graphically  representing  equation  (A. 2)  results  in  Figure  A2. 


0  t  -  €  t 

Figure  A2.  Step  Time  Delay  Function  J(t  -  £) 


Substituting  into  the  general  time-delay  integral  the  particular  step 
function  in  terms  of  the  non-dimensional  variable  £ 


w(t)  = 


ft  -  t  dm. 


ft  dtu 

j  1ard£- 


Therefore,  simplifying  equation  (A. 3) 

u»(t)  ■  wQ(t)  -  wQ(t  -  t) 


(A. 4) 


(A. 5) 


where  u>Q(t)  is  a  generalized  function  of  time  and  wQ(t  -  t)  is 
functional  time  delay. 


APPENDIX  B 

RUNGE-KUTTA  PROGRAM  OF  THE  MODAL 

s' 

AMPLITUDE  WAVE  EQUATIONS 
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GARY  ^  1  COON  A  LO  SYSTEM  OF  QPOINAAY  SECOND  3P0ER 
DIFFERENTIAL  EQUATIONS  FOR  NON-LINEAk 

combustion  Instability 

GCYEFMNG  EIUATIQNS-CASE  NO  TIME  OELAY  (J*Q) 

tl<  T»Fl»F2*(»i»GE*Fl®»F2P#Gl,?#G2r*  )»F1P 
S2(T,Fi,F2,<il»G2»FlP,F2P,GlP,G2Pi1!F2P 
£  T(  T  ,F  l  ,F  2,GI»G2»F1P,F2P»GI<*  »G2°) »GIP 
E«(T  #F1»F2#GI, G2#F1P*  F2P#G1P#G2P  >*G?P 
£M(T  ,Fl»F2»GWG2»FiP»F2P»Gla.G2P>*  (-F 1- «9 Afi *F IP«2 . *£ P 
1  •1*(F2*F1P  *F\*F?.P*  G2*31P*G1  *3 2? > -2. *N*£P* M9 A F * ( F 1*F 2 
2*C2* il))/( 1. ♦**£?> 

£o<  T,Fl.F2»Gt,G2»Fl*».F?P,Gl»  ,32P  >»(-4.*F  2-*3  AR*F2PTP 
1  *  1*(G1*GIPFI*FI?  )-0.S*':P*«48A"*N*(G1  •*  2.-F  1*»2.  )  >/(  1. 
2*A.*N«E°) 

£7<  T,Fl»e  2»Gi»  G2,FIP,  F2P.3I*  ,G??)s<-3l-*3A**3lP-2.*£P 
l*l*C  12*FIP«FI* G2P-Gl*F2P-F2*.UP)-2.*N*t>*  WBA^*(F1*G2 
2-F2*GD)/(  l.«X*£P) 

S2<  T#F1»F  2»Gi.G2,FtP*F2*»,GlP,G2P>*(-4.*G2-w34R*G2P<EP 
i  *I*(  ll*Fl?«Fl*GtP)*.HA:>*N*£P*(Fl*Gl>> /(  l  .♦4.*X*£P  ) 

PEAL  K.N.I 
L  8 1 
T  8  J. 

5-:u  initial  conditions  and  constants 

E^-OAOE"'  TERM  i3$IL0N 
WNAP-STEAQY  STATE  burning  rate 
win? action  inoex 

ia-COFRSCT  ION  TE«M  (BAFFLES  NOZZLES  »E  TC.) 

I-(JAS  0*NAMIC  INOEX 

Ft, F2. 31,  G2-M00AL  AMPLITUDES 

A£A3  (5,W>  £?  ,  W34Rt,N,X,l,\MA  X 
10  FCRMAT  (  5C  l  A . 4  »  I  i  ) 

REA  0  <0»t?>  Ft,  F1P,F2.F*P,G1»GIP,G2,32» 

Is  r JP  N  A  T  C3F|0.4> 


>•  I  A  D  STEP  SIZE 


READ  (5*1A)  N 
14  FCr^AT  (Fl^.S) 

XRITE  ( &»?0>  £P,NSaP.N,K,1,H 
2 C  FO'NAT  C1X.4X, ‘EPSILON*. HO. 4,/, IX, 4X, ‘STEADY-STATE 
1  BURNING  PAT£(  N  AAR  )*, FI  0*A#/»IX»4X,‘ INTERACTION  INDEX 
?-N»  ,F  10  .4.  /  .lx,  4  X.  ‘CORRECTION  YA* 1  A3 1  £-  K  •  .  F  1 0 . 4  ,/  ,  l  X 

3,  AX, ‘GAS  0  Y  N  AMI  C  N0N-LINEA»  I NDEX- (  I  >  • .  FI  0.  4, // ,  l  X,  4X 

4,  ‘STEP  SIZE*  ,Fi  ?.<♦»✓/> 

m  F I  T  E  (  G,  J  0  )  Fl,FIP»F2,F2P,31  ,G1P,32»32P 


o  o  o 
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3C  FCi^AT  CIX»'I«1MU.  C3N0I  M  0^  S*  »//»lX»4X»  *F  l  *  »F  H.4*/ 
l*lX*!»X.*ari/0T».F10.4,/*U,4X#»f2»  .FI  0.4#/#  IX.  4  X»  »0F2 
2/:T»  » FI 0.4  */*lX,4X,  *31 *#Fi0.4#/#lX.4<»  *3G  l/3r »»F13.*'/ 
3  *  IX  ,  4  X  »  "32  »  »F  l  .  4  ,  /  , l  X  ,  4  X  *  •  3G  2/0T  *  •  F 1  i  .  4  ,  //  ) 

ARTE  (6#40> 

4a  FCRMAT  ( IX  »  3X »  *TI'*E*»?X*  •  FI  * »  1  OX*  *F2*  »  1  OX  »  '  G  l  •  *  9X  # • G?  * 
1  *  1 5  X#  *  AC  J1  •  #!  JX,  *  ACP2  *#  ?X,  'ACP3»*  9X#  *  AC**  '»/> 

PUNGE-XUTTA  A  L5QM TH* 

3C  IJ3  J  *  l.NMAX 

c 

°l*H*El(T»  FI,**  2,  <11  ,  G2»  F1P,  F2-»,3  1P,J2>  > 

31*-*  *£2(  T»F1  ,f  2  »  31  ,  G2  ,  FI  P»  F2P  ,  Gl  P,  i?»  ) 
•U*uUU<T,Fl,F2>3l,'j2,FlP,F2P,Gi»#3  2®J 
S  Is*  *S4(T*  F  t  »F2>  31  ,G?,FlP,F2-»  ,U1P,32?  1 
Tl*rt*t5(T*Fl*F2»jl#G?#FlP#F2P#GiP#.l?») 
UI»H*E6<T»F1»F?»3I»G2*F1P#F2?»510»>j29> 

/ 1  *  H  *  Z 7  {  T,  F 1  ,  F  ? ,  5 1  ,  32  ,  F  I  P  ,  F  Z?  ,  3 1 P  ,3  2P  ) 

4i*M*£J(T,Fl»F 2»GI,G2*FIP,F2?  ,G1P,32°) 

C 

?2  »H*  tl(T*H/2..Fl4Pi/?.,c24  qi/2.  ,Gl*»  t/2.#G2*S  1/2.»F1P 
l«Tl/2.»F*P4Ul/’.*Gl»4«l  /  2., 32 ?4*  1/2.) 

Q2*H*£’(T4  H/2.»  FltPl/2., F2«21 /2.#Gl 4fi 1/2., G2«  s  1/2.»FIP 

!♦  T  ;✓  2.  #F?P«  U 1/ ?.*  Gl  »»♦  xi /2.,3>p  «4i/ 2.  ) 

■J2=H*£3(T*  H/2..Ft4Pi/?.,F2.  jt  /2..GUR  i/2.,r>2*si/?.,rlp 
!♦  T1  /?.#F2P  «G1/  2.  »G  IP  «  y  1/2.  ,G?P4  « 1/2.) 

i2*H*£4C  T4  M/2 .»  Fi 4P1/2., F’431  /  2. *  G1  1/2. » 12*  SI /?. *F1P 
l4Tl/?.*F?P«Ul/2.»GlP4tfl/2.»3  2p«<it/?.) 

T«  =  J4E5(T«H/.'’.#Fl4?1/2.,F24Qt/2.,Gt4Rl/2.*j2  45i/2.#FlP 
14  Tl/2  .,F;P4'JI/2.,G1P*  VI/ 2.#  32P4W1/2.) 
j2*H*FSCr4H/2.#Fl4Pi/2.,F24Si/2.,Gl*'«l/2.*G24Sl/2.4FlP 
l4Tl/2.»F'?4Ul/2.,GlP4Vt/2.,r,2p4Wl/_>.) 

/2*M*EF(T4H/2.,Fl4pi/2.,F24  31/2.,Gl4^l/2.,G2*Si/2.»FlP 

14TI/2  .  ,F??4*W1/2.,GIP4V1/2.,G2P4  wl/2.  ) 

A2*H*  t?(T4  H/2.,Fl4  Pl/2.,F2«  <n  /  >.*Gl  **  1/2.  #32*  S1/2.»FIP 
l4ri/2.>F2a4Ul/?.#31P4yi/2.,G2p4Wi/2.) 

% 

w 

°I*H*cl(r4H/2..  Fl4P2/2.,F?4  !3!/2.  .QUs  2/ 2.,  3  24  3  2/ 2  ••F  IP 
1 4  J2/2.,F2?4U2/?.*atP4V2/2,#G2P4W>/P,  ) 

;S*H*E2CT4H/2.  »  Fl  4P2/2.,F24(32/2.,Gl  4«  2/2.#324  52/2,#F  IP 
l*  T2/2.»F2P»U2/?.4GIP4  y2/2.,G2P«*w2/2.  > 

J(  T#  H/2,.  F 14  P2 /2.,F  202 /2.  »  UUP  2/ 2.,  G24i2/2.,F  IP 
l«I2/!..F?p4{j2/2.,Glp4/2/2.,G2P4W2/2.) 

i  3*rj  *  £  A(  T4  H/2..F14  P .?  /  2  • » F  2*  32/2. ,Gl 4« 2/ 2 . , G2«  32/2.  ,?IP 
1  4  T2/2.  ,F?p4U>/  2,  ,GIP4V2/2.»G2P4  W2/2.) 

T3»H*iSC  T4  h/2.^  FI  4»  2/2  .,F  24  32 /2 . » Gl 4  R 2/2 . , G? ♦  >2/2.,FlP 
i*  I* /%*F?P4L2/?.,G1P4/Z/Z..G2P «*2/2.  ) 

U3*4*Ej(T*H/?.  ,  Fl4P2/2..F24Q?/2.#Gl4  *  2/2.  #32*  3  2/2.  »F  t«* 
l*T2/2.*F?P4U2/  !.,3lP4y2/2.4u2P4W2/2.) 
yj*H*EF(T.  H/2..  F14P2/2.,  F24  32/2.,Gl4-?2/2.  ,q?4  S2/2.,FIP 


-  A-**  •* 


1  ♦  TO  /!  .  *F2?  *U2/2  .  *  G1P*  V2/  2.  *32?*  W2/?.  ) 
.*3=“t*E3(T«H/2.*fi«P2/2.*F202/2.*Gi4P2/2.*G2*S2/2.»FlP 
1*  T?/?.oF2°«U2/?.»GIP4V2/2.#G2P*  W2/2.) 

%=*4*£l(T4H»Fl«P3#F2*j3#Gi4SS»32*s3#Flp*T3#r2:,«u5*GlP 
!♦ *3*  32?4*3) 

S4a  .-)*£2(  T*  H#F1  4P3#  F2  42  3#  Gt  <P  3/G24S  3#  F  l»4  T3#F2P4  US»G1P 
l*  V3  *32P*  W3  ) 

MA*H*£3(T4H»Fl4P3*F24.J3,3l4h3»G2*S3»rlP*T3*F2p4U3»Gl:» 

l ♦ 3  *  J2P4  W3 ) 

34*H*£<*(T«H#FMP3»F2  433»Gl4P3»G24Sl#Fl»4T3»F2P4U3»Gia 
1  *  w  3  *  G2?  ♦  VJ3  )  „  .  . 

TA=H*£5(T4H»Fl4P3*F24  3  3#Gl4P3#G2-*S3»rl?«r3#F2?4U3»GlP 

1  ♦  k3*32P4‘<3) 

U4  =  H*£GCT-*H#pl4P3»F2  4  3  3»Gl4P3»‘>24S3»FlP4T3»F2P4U3»Gla 

t  ♦ i\ *3  2°« ^3  ) 

/4-H*-:?  <  T  ♦  H*  Fl  *P  3*  f  3*3  1«*  3»G2«S3*FIP4  T3»F2?4L3»GIp 
i*  il'  32P«W3  ) 

rf4sM*£e  <  T*  H  *Fl  4  P  J*F2«  3  3*Gl*P  3*G  2«S  3»"  lp*  T  3*F  2a*  U  3*G1P 
1  ♦  v  3*  3tp*  #3  ) 

Fl=Fl4((3i42.*P242.4;*34P4)/6.  ) 

r  g  s  -  ?  «  ( C  3 1  «  2 .  *32*2.  *5  3*  3  4  )  /  o .  > 
Gla|ji4((Pi'*2.*P242.*<?J4K4  )  /  3.  ) 
i2  =  j2  4((5l4?.*3  24  2-*534ai*)/o«> 

Fl°=Flp4<( ri42.*T2*2.*T34 T4) /5  •  > 

F2r,sr2P4<C'Jl42.*U242.4U3*U4)/6.) 

5  l p  =  i  l  P  ♦  <  <  X 1 « 2  .  *  ^2  ♦  2 .  *  3  ♦  M  > /  o  .  > 

3eJ=G2p4((  WI*2.*W242»*«<3*  A  4 )  /  j  •  ) 

}£-»i--Pl,3«E;PS*(-Fl*F2*3i*j24  3.5*(“»lp*j2P4FlP*F2P>) 

ACp2="F2p«  £  =* S*- <  - G.25*(  0  i  *31 -r  1 *F l )♦ 0.  2?*(F1P*F1°*GIP 
l *ci p  > ) 

ACP3-“'jlp4E0S*("Fl*32»c2*Gl4  3»’*(FtP*j2P“GlP*F2P)) 

4Cp4s-32p4£pS*f0.5*Fl*3i40.5*clp*GlP) 


T  *  T  ♦  H 
L  *L  *  l 

IF  (L  .E3.  3)  3Q  T3  110 
30  TO  103 

110  *P1T£  (3*50)  T  »Fi,F2*31*32*A:Pl*ACP2»4Cp3*A':P4 

30  FCSH4T  C1X  *F3. 4  *  JX*  FI  *.  5*  2X*Fl  0.  5*  3X,  FU  .5*2X*F  10 .5»3X 
l*F10. 3* 2X»F 10.3*  3X*C 10.5*  2X*F 10.3 > 

Lsl 

103  CONTINUE 
CALL  EXIT 
tNO 


iwiuf  m  wii*i  w  pi,  . .  w  ^  >4 ... . #fT,Tvrnwv^TTTT!yw..  w 

'  ■  -  ,  ^  ^  ^  ^  t  *  i  ^ 
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APPENDIX  C 

RUNGE-KUTTA  PROGRAM  OF  THE 
PERTURBATION  EQUATIONS 
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G  A*  V  3  MCDONALD  M0N-LINT4A  C0M3USTI0N  • -3U  AT  I  0N$  DERWEO 
FROM  THE 

TWJ-*4H!48LE  PESTUfiflATlON  ME  ThDQ 

GOVERNING  EaUATIJNS-CASE  -  NO  TIME  DELAY  <J*3) 

£HTj  U#3l  *42*32*4  3#3  3*44*94)*<-0.5*iI.33*Al-0.5*K*8l 
l-C.5*l*(4l  *4  3*91*93<4?*A44  92*84)-0.5*N*rf3AP*<  3I*43-Al 
2*E3*  32*44-42*34  >)*EPS 

E2(T*  41*31 * 42*3  2*4  3*3  3*44*34) *(-0.5*5 IG3*  8l*  0.  5**<  *A  l 
1-3. 5*1  *(41*933 1*43*44*32  *42*3  4  )-J.5*N*N34f»*(Al*A3*31 
2*S3*42*-A4«  32*34  )  >*£P$ 

E  St  T**i*3ljA2»  32*4  3*9 3*  4  4 *3  4)  *(  -  0.  5*5  JG3* 4  2*’ 0*  5  **  *3 2 
l -0.5*1  *<4  1  *A4<  ?  1*3  4-42*4  3-32*  3  3  >-0 . 5*  N*R3  Al*  *(  44 *3 l -At 
2«£4M2**3-A3*32>)*E(,S- 

£4(r*4l#3t#  42*  3  2*4  3*j3#A4*34)  *l-0.5*SIG3*82«  0.5*K*a2 

1- C.  5  *1*04  *4  144  *31  4*3*32-42* 3 3  >- ) .5  *M *R3  A R  *(  4  1  *  A4*  8 1 
2*64-42*43-32*33  >  >*  £PS 

E5<  T*  Al  #31  *  4,!*  3  2*4  3*93*44*3  4)  *(-0.5*5  Iii9* A  3-4.  *K*J  3 
'  1 ’0.125*1  *<  42*A2-32*,32*3i*3t-4  1*4t)-0.  12 5*W34R*N*(  41*31 

2- A2*  32  )  )*ZP5 

EE(T#At*31 *42#  ?2#A 3*33*44*94) *(  -  0.  5*5  ID  3  *8 3*  4.  **  *4  3 
1-C.  2  5*1  *(A  2*  32  4  1  *3  1  )-T  .06  25M  34R *N* ( 4 2* 42 -92 *32-  41  * A1 
2*  El  *  31 )  )*E  PS 

EF(T*4l>31»42*32»43*33*A4*34)3(-0.5*5IG3*A4-4.*K*3  4 
1-C.25*!  *(3  1*324  1  *42)-0.125*W94s*m*C  At  *32  *  42  *31  )  )*EPS 
t  3  ( T  *  4  1  *  3 1  *  42  *  32  #  A  3*  3  3*  44*  3  4 )  * (*  0. 5 *51  33  *9  4*  4  ,  *K  *4  4 
1*0.25  *1  *CJ2*4U  4  2*91  )*0.125*<  8AP*N*  (41*A2-31*3  2))*EPS 

D  MEN  SION  Ft  C  5  »>  0  )*  F2C  5  30  )*  G1  (  500)*-3  2(  530) 

3  •  A  •.  X  *  n  •  1 
L*1 
r*o. 

BEAD  INITIAL  CONDITIONS  AND  CONSTANTS 

EPS-OROE-?  TERM  (EPSILON) 

WNAh-SfEA  DY  STATE  BURNING  RAI¬ 
N-INTERACTION  INDEX 

K-C3RRECT  ION  VAHIA3LE  (9AFFLES#  WALL  LI  NI  NGS  '.NOZZLES 
» £T  Z  *  ) 

1-34  3  DYNAMIC  INDEX 

Al*31*A2*32*A3*d3*  44*34-M004L  AMPLITUDES  COEFFICIENTS 

READ  <5*10  )  wJ4R#N#K.NM4X*EP5*  I 
1C  F0PM4I  (  3F10.4*  I  J*2F1).4  ) 

READ  C S *  1 P  )  41  *  3  1*  42*  32*  4  j*3i* 44*34 
12  FCPM4  T  {  ?r  1  0.4  ) 

C 

C  FEAO  STEP  SIZE 


non 
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*£0  (5,14)  H 
14  FOMAT  (F10.5) 

Sni*V3AA/£R$ 

rtPlTi  (6,20)  $IGd,rf3Aft,N#K»H»EP$»I 
20  FCPMAT  (IX  ,4X,»  SIGMA  8  \r  •  #  F10 . 4  #  /#  l  X#  4  X#  *  ST£  AD  Y  $TAT£ 
l  3U4NING  »AT€-(M9AF)*»f  10 .4 */ , tX t4X t « I  NTS *ACT I OM  iNDtX 
2-(M  )•  ,F10.  4,/,lX #4X,'  C0» ^ £C  T  I  QN  VAR  I  ML E- (  X)  * #F10  . 4  ,/ 
S, IX  #4  X,  'STEP  S12£-(H)  »,Fi:.4,  /,1X,4X*  •  C  P  S 1 1. 0  M  “  (  E  P  5  )• 
4,F10.4*/»IX*4X#'GAS  DYNAMIC  IN D£«- ( I ) • *  FI  0. 4 *  // ) 
rfftITE  (6*30 )  Al, 31, *2,32,43,33,44,34 
so  format  (i<«' initial  C3N0tT!0NS',//,ix,4x,  »ai*  »Fn.4,?x 

l»«Bl»  ,FIQ.4>/,  IX,4X,*  A2',F10.4, 3x,  »32*»F10.4,/,lX,4x 
2#  'AS*  *F13«4,8X,  •  83',Fl0,4,/,lX,4X,,A4,,Fl0.4,'iX»'34* 
3»Fl J.4,//) 

*  s  I  T  £  (5,40) 

40  FCMW  (  1>:  ,3X*  •  T  IMS',  7X*  »Ai ',  >X,  ♦&!',  12X,  *A2' »  IIX, '32  • 
t*i:X,'A3'»llX,'33*»lOX,'A4»,lCx,'d4'*/> 

KLNGt  <UTTA  A  LCDn I  Tr|M 


oc  no  j  *  i,\*ax 

Pl*H*El(T,.4l,ai,  A2,3>,  A  3*  3  5,  U,34> 
0l*H*E2CT,Al,3l,A2»a2,43,33*A4*$4> 
ft  1  »H*-  £  3 <  T,  A!  #31#  42,32#  A  3, 8  £ ,  4  4 , 34  > 
5l*H*£4(T,Al,3l,A2,ci2,4J,3  5,4  4,  34) 

Tl*H*S5(T,  Ai,  31,  42,32*  4  3#  33, 4 4, *4  ) 
Ul*H*E5(T,4l*il,42,32,AJ»,3S,44,34) 
Vla'N*EF(T#Ai*3t*42,8?»A3»'3S,A4,34) 
a l* M* £3 (T, 41,31, A2, 32, 43,33,44,24) 

P£*H*  £l(  Tt  H/2.,  41t?i  /2.,3tt  SI  /2«,  A  *4*  1/2.  »s|  2*  5l/2.  ,*3 
ltTl/2«»3  3«Lt/2*»4  4*Vl/2,»J44i(l/2.) 

C2»H*£2(T*H/2.»  Al*Pl/.!«'4i*3l  /2.,  A  2  *4  l  / 2  3  2  *  it/2.,4  3 

It  TI/2.M3*  Ul/2.  »A4*  Vl/2.#34  4Wl/2.  ) 

4  2  =  N*S3(  T«  i-l/2.,  A14P1/2  .,3  14  01  /2.,A  .'tR  l/2.,3  2t  Sl/2.,4  3 
1*  T1  n.'il*  Ul/i.,  A4t  Vl/2.  #  34  4  ^  l/>.  ) 

S2=H*£4(T4M/2.»  AI tp 1/2. *31*41/2.' \?4R1/2.,R2«  51/2.# A3 
l4fl/!.,35tUl/2.,  A  4, V l / 2. *3  44  U/2.) 

r2*4*£5<T*H/2.,Al4?l/2.,Bl4(H/2.,42*5  l/2.,32*  Sl/2.,AS 
It  TI/2.»8  3tUl/2«  »  A4tVt/2.#.3  4«w  1/2.) 

U2*  i*£GCT*H/2.»  AitPl/2.*3l4Ql  /2.»  A>4i?l/2.#  3  2*  51  /2,,A3 
l  4  Tl  />.,3{t  Ul/2.  ,44* VI  /2.,34<W  1/2.  > 

V2M*£  ?<  T  t  H/?.,  A  1 4  (»l/  >.*11*  0l/2.#42«P  1/2.#  324  51  /2  .  *  A3 
1'  II  /2  .  ,3  3t  Ul/2.  ,  A4t  vt  /2. ,  34  ♦  *  1/  2.  ) 

N<sH*£d(  f*  H/2.»  Altpi  /?  .,3lt  Qt  /2.  ,A2tP  l/2.»  3>4  >1/2.,  AS 
14  Tl/2.,3  3'  Ul/2.  ,  A  44  V 1/ 2.  #84  4  W  1/2.  ) 

PS«4»clcr*H/2.,  Al4P2/2.,31tQ2/2.,A24R2/2.,3  2t  $2/2.*A3 
It  T2/*.,3  3«  U2/2. , \44V2/2.,d44W2/2.) 

Si  =  H*£2( T«  H/2. ,  Al 4p2/ 2. ,  31  4 02/ 2. , A2tB 2/2 . ,3 2' 5 2/2 . * A  3 
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t?/>.*33«U2/2.>  A4«V2/>.,  84  *142/2.  I 
*  Ja4fr£3(T«H/C.»  Al«P2/2.»3l«a4/2«»A2«*2/2«»324  $2/?.* A  3 

14T2/2..8  34  U2/2..  4  44  72/2.,  9  44*2/2.) 

S3*H*E4<T4H/2..Al4P2/>.,9t402/2.»A2492/2.,624$2/2.,AS 
^T*/?.^*'^/?..  44  4  72/ 2.  ,94  4  *2/2.) 

r3*H*t‘5(T4H/2.»  4l«P2/2.,dl*  32/ 2. »  A  24  *  2/ 2. ,  9  24  *  2/2 .  ,A  3 
14 72/2. .3 34 J2/2. .A44V2/2..344W2/2.) 

U3»H*t6(T4  H/2.  »  At  *P 2/ 2.#  3  14^?/ 2., 4  2* R 2/ 2.  .824  S 2/2.. A  3 
14  I2/2  ..334  U2/2.  .  A4  4  72  /2.»84*42/2.  ) 

V3*M*Em*  H/2*.  A14p2/2.#3l4Q2/2.,A24R2/2..32«s2/2..A3 

14T2/2..334 j?/2..A44V2/2«*944M2/2.) 

W3*H*£3CT4  H/2.«  Ai4P2/2,,8l4Q2/2.»  A24R2/2..  3  2a  S2/2.»A3 
14  T2/2  ..3  34  U2/2.  .A44  72/ 2. #04  2/2. ) 

C 

P4«H4E1CT4  H.A14PJ.  8i4Q5.  A24p  J,8  2.$S.A34  T3.3  3*  U3.A44/3 

1.644  .13) 

34*H*i2(T4H.  At4P3.3l4«3»  A24f-3»32*S  3.4  34T  3.334U  3.A44V3 
1.  64. A J  ) 

P4»H*£3  tT4H.  A14PJ.31  4Q3,  A24P3.924S3.A  34T3.3  34U3.A44  73 

1.644  *5 ) 

S4«H*£4C  Tt  H.  A14P  3.3  14  03  .  A  2*  9  3. 92  +  5  3.  A  3*  T  3 , 33  4  U  3,  44  4  7  3 
t. 64*^3) 

TA«H*£5(T4H.Al4P3.8t403.  A24R3 .324S3.  A34T3.3  34  U3.A44  7  3 
1 , 34 4  W3  ) 

«4=H*£,i{T4H»At*/*3,9l4j3»4.24fi3.924S3»A34T3.i34^3.A4473 
1.  £4  4  14  3  > 

74«H*£7(T* H,  41  4  P  3.  3140  3.  4  2*  4  3.6  24$  3.4  3*  T  3. 3  3.  U  3  »4  4*  73 
1 . 84 4  *  3  ) 

44« H *£3  (  T4  H.A1  4 P  3.  81  43  3.  A 24R  J,  3  24  5  3.  A  34  T 3.1  34  US  >4  44  7  3 
l.d<*.  43) 

C 

Al»4l4t(Pl«2.*P2A2.#p  J«p*  >/$.  > 

3  t« 114 (<  31 4 2.* 02 42. *0  3 434  )/ 6.  ) 

A?,  *424(01  42«*i.  242«*934f;  4  )  /  6  •  ) 

92‘624(<it4  2.*S2«2.*S34 S4)/o.  ) 

A3*434(<n*2.*T24  2.*T34T4)/6.  ) 
33«1l4t(Ut42.*U?42.*Ul4U4)/5.  ) 

4  4*A44((7l4  2.*7242.*73*74  )/6.  ) 

3  4*  144 C( 71 4 2^*4  24  2, *7 34 44 )/6.  ) 

C 

T»T4  ti 

C 

Fl(J)«At*CQS(r>43l*5I^(T) 

Git J)*A2*C0S(T )4J2*S!XCT  ) 

F2(J )«A3*C0S<2. *T)*93*SINC2.*  I) 

G2(  J)«A4*:0$(2.*T)4d4*SI  M(2.*  T) 

C 

L*t4  1 

IF  (L  .£3.  3)  10  TO  110 

10  TO  100 

110  XFITE  (6. SC)  T.  Al.6l.A2.92.A3.63.A4.94 


50  FC**U  (iX*F8.4»  3x»FI0.5*2X.Ft9.5»JX.F  13.5.2X#FlQ.5*SX 
1»F13.5.2X.F10.5#3X*F10.5.2X.FIQ.5> 

L»l 

iOC  CCMTINLE 

1*1 

t«0 

rtfl  TE  (5**1) 

41  F3MU  (1X.///.1X.3X*  »TI4£»*>X#»Fl#aOX#»F2»*llX»»Gi' 
1#  12X.  *G2».  /) 

oc  2 o c  j  «  l.^yxx 
T«r  ♦  h 

L*L«  1 

If  <1  .EO.  3)  GO  TO  210 
3C  TO  230 

210  IT2  (5*300)  T.F1  ( J).  F2 ( J) . U < J >.32<  J) 

500  Fa*.  4  \  T  <lX.F2.W3X,Fl0.5*2X,.-i0.5#3X.F10.5»2X»Fl0.5) 

L  *1 

20  0  :c*jtinle 

CALL  EXIT 
ENO 


i 


APPENDIX  D 

PROGRAM  FOR  EXACT  SOLUTION  OF 
STANDING  WAVE  CASE 
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C  G  if  H  >1  COON  A  LO 

C  ANALYTIC  SOLO  T ION*  STANDING  WAVE  CASE 
JHENSION  FU  5  20), F2C  500) 

PI  =  3.  1415D2S 
H*  0  •  1 
X„*0  •  0 
NMAx»500 
E  F  3  *  J  .  I 
W9AR*C.l 
Lal 
C 

C  READ  IN  GAS  DYNAMIC  INDEX  AND  INTERACTION  INDEX 
C 

READ  <  5#  5 )  J,XN 
5  FORMAT  CI5*F1C«5) 

C 

C  INITIAL  CONDI TIONS 
C 

4 E AD  (5,1)  >  C 1 »  C  3,  PHI  1 »PHI 3 
1C  FORMAT  (4F10.4) 

C 

WnlTE  (6,15)  I 

15  FORMAT  ( l  X » 'GAS  DYNAMIC  !NQEx*,l5> 

WRITE  (6.15)  XN 

16  FORMAT  (IX, ’INTERACTION  INDEX » *?1 . Z >! // ) 

WRITE  (S»20)  Cl #C3#PMI  l,PHI 3 

23  FORMAT  (1X»»  INITIAL  COMO  I  TION  5  //.  1  X,  4X,  »C  I  '  *F  ID  .  4  ,  t 

1*  LX,4X,*D3*  #F  ID  «4»/,lX»4X'  *  PM  1 1  •  ,F10.4*/,  lX*4X»,e*HI3* 
2#FlO*4,/> 

WFITE  (6, 30) 

30  FORMAT  (  IX  »2X*  •  TIME**  L4X**  Al*  »13X»  *  A3*»13X»<31*  ,l2x' 
i>*33*  ,  15X,  »Cl*#1ZX,*C3  *#!1X,*PHI1*,19X»  'PHI  ?*»/) 

U 

JO  l’DQ  J  a  1  ,  N ‘1  A  X 
X*=X J^t^S 

IF  CX,\  .£3.  D..AN0.  1  .E3.  1)  GO  TO  f>o 
IF  (  XN  .ED.  AO.  .AND.  I  •  E 3 •  0)  GO  TO  50 


C  \C  GAS  DYNAMICS 
C 

50  CONTINUE 

fi«  =  l  .0 

S*<;. 35353 •XN*EPS*(-l.)**J) 
T*I.-:x°C-Q.5*w8ARwxJ} 

UO*T 

V*l. /CQS(U) 
w*EXP("0.5*w3AR*XJ) 

C  1  s  M  *  V 

Z*(SIN(U)/CDS(U)  ) 

C3  =  (W/2.’3234)*>Z 
Phil *PHI1 
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>*Pi 

jj  n  2oo 

* 

M 

C  NO  CONiUSTllM 
C 

*0  CCNTINUE 
p  «  *  i  • 

1S« <£P3*(«1,  )**J)/(?.*43  4S)  • 

CT*i»"SX?(*0«5^w3AR*xj) 

ci*es*cr 

Ci/«2, /<rxp  (CUM  l,/CXPCCU)) 

C***:XP(-;?,  5*»HMu) 

Cl *CW*CV 

Phll*PHll 

PHI3=.2«*»HJt-;j*(»j 

2oc  continue 

A l*Cl*CC$<  PMli  ) 
ii=ci*: iNtPHn  > 
am:  ucj$(?hj  j> 

•3  J  =  C  J*SIN(  PHJ J  ) 

F1(J)  =  AI*CJS(XJ)*31*$Iv(XJ) 

FZ<  J)  *  AJ*C0S<2.  *  XJM3  5*SI<\(  2.  *XJ  ) 

IF  (L  «EJ.  ?)  10  TO  5'iv 
c=L  ♦  l 

u  rj  iu.) 

500  TM-:  =  XJ 

,  -iXcrj'-'JX'  li>*5'3X'Fl0*0*tx'  F»0.5»;x.rn.5,5x 
lIj1  5°X#  °*3#4X'f’l0,5#4X'r  1  )#  5>  * 

no  CCMUnUE 


X  «  =  3  • 

L  =1 

^f*IT£  (MaOOO) 

^6100  FOPNAT  UX, ///.IX,  IX,  •?!*£. .l(ix,«Fl»,15<,*f  »»,/, 

JC  J 0  0 j  *  1,‘jmAX 
X„*X J«H 

if  cl  .eo.  2)  ao  ra  ro-i 

L  *L  *  l 
30  TO  300 
7 00  THS-XJ 

_.A1.  **!  (6^>000)  T1M£,F1<  J),F2<  J) 

F  000  FQMAT  <  1 A  *FF.4  #8X*F1').5  »F  X»- 1  0.  5,/ ) 

30C  CCNT I Nw£ 


igMIij mjk 


r  i  o 


all  -:xit 
ho 


.Uifc.  fc.  :w.U  fafr  fe'  | 
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C  3AF 1  H  MCOOMALB 

C  ANALYTIC  S3LL  T 13N* TH A/ELL1NG  mAVE  CASE 

DIMENSION  F1C5)0)*F2(  500>*U(  SCO)*  i  2(  50  0) 

P  Is  5  •  1 4 1  *a 9  26 
H*3  •  l 
X. »0  *0 
NMVX*500 
EF5  =  0*  1 
A3  AS *0  •  1 
t*l 
C 

C  PEAO  IM  0  AS  DYNAMIC  INOEX  AMO  INTERACT  ION  INDEX  , 

C 

*£A0  (3*5)  I  *  X  N 
5  FO-'-IAT  (I5*F10»5) 

C 

c  initial  conditions 
c 

flcAO  (5*10  )  Ai*  A  2*  A  3*  A4 
IC  FCM  AT  (  4F  10. *4  ) 

AfilTZ  (6*15)  I 

15  FCSM  A  I  (IX* 'GAS  DYNAMIC  I NOEX • * I  5 ) 
hFITE  (6*16  )  X M 

16  FORMAT  (IX*  ’INTERACTION  I NOEX  '  *F  T  .  Z  */ It ) 
nnllE  (6*20)  Al*A2*Aj*A4 

20  FORMAT  ( 1<  * ' I N i T  I  A  L  CON 0  I T ION  5 * * // *  1 X  * 4X *  •  A  1 • *F 10 . 4  * / 
l.lX»4X**A2<*F10.4*/*iX»4X*‘A3**FlC.4*/*lX*4X*tA4* 

2*  F l ). 4*/) 

(S.iD) 

30  F  C  f*  M  A  T  ( IX *4X* *T  !*£«* 14X *  * Al*  * 15X* * A2 » *  15 X* *  A  3* *  1 TX 
1**44**/) 

DO  no  J  =  l*M-4X 

Xja^j4£P3 

IF  (XN  0 • »  A  NO  •  I  .£3*  1)  GO  TO  60 

IF  (XN  .£0.  20..  AND.  I  .£3.  ))  GO  TO  5  0 
C 

C  MO  GAS  DYNAMICS 
C 

50  CCNTINLE 

S«O.FO?ll*tPS*XM 

T*1.-EXB(-D.5*W3AR*XJ) 

U  =  S*  f 

V*t./COSCJ) 

A*Z X»(-0.5*W3AA#XJ) 

A2*«*  V 

/*S1  UU)/COS(U) 

A4*( A/l. 62642) *Z 
GO  TO  200 
C 
G 
C 


NC  CO-HUSTION 


60  CCNT 1 MUE 
CS*E»»>/W34f» 

CT*l  .-£XP(  -S.5*h9AR*X  J) 

Cll*C  5  *CT 

C^2./<EX?(CU)«  l./EXP(C'J)) 

C***X?<«0.  5*w34«*XJ> 

A1*C'W*CV 

C2*(  EXP(C0)-1./EXP(CU)  >/(tXP(  CU  >«t./£XP(CU>  > 
A3*(C*/2.)*C  2 
3C  TO  200 
200  3i*-A2 
32*41 
33*  *44 
34*41 

Fit  JJ*A1*:0S(XJ)431*S1M<  XJ) 
ai(J)*A2*C0SCXJM3  2*SIW(Xj) 
F2(J)*A3*C0S(2«*XJ)H3*$I.M(2.*XJ) 
32(J)*A4*C0S(2.*XJ)4  0  4*S1M(2»*XJ) 

IF  (L  .E3.  2 )  30  TO  300 
L*L'»  l 

go  n  loo 

500  TH£*XJ 

/•FITE  (6»  5000)  TI*t#Al»  42#A3#  A4 
5000  FCFM4T  (IX  #F7.4  #7X#  (4  (Fi2.  6*3  X  ) )  ) 

L*1 

100  CONTINUE 


X**).  C 
L  *1 

aPlTE  (S#6000) 

6000  FO^AT  (1  X# ///#  IX#  IX#  »TI  *£  •  #t  2X#  *Fl  •*  15X  *  »F  2' *  16X  »  »3i  * 
1# 14  X  # *32' # /  ) 

O’J  30  C  J  *  1  # 'NH  4  X 
XJ*XJ«H 

IF  (L  .£3.  2)  30  TO  7  00 
l*L«  1 
3C  TO  300 
700  T  I M  2  *  X  J 

WPIT^:  (6/7  000)  TlM£#Fl(j)#F2(j)#0l(J)4S2(J) 

7J0/)  FOiHA  T  (IX  #F7.4#3X#F10.5#?X#“  10 .5#8 X# F  13  .5  #7X  #F  13 . 5  > 

L  *  l 

3oc  continue 
call  exn 

£F0 
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PRESENTATION  OF  ACOUSTIC 


PRESSURE  CALCULATIONS 
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ACOUSTIC  PRESSURE  DERIVATION 

To  calculate  expressions  for  acoustic  pressure,  recall  equation 
(2.48)  which  stated 


p  =  P  = 


ft  *  ^  e2^2  +  2u  H*  +  V$  •  ^ 


(F.l) 


This  equation  represents  the  unsteady  state  deviations  of  acoustic 
pressure.  When  expanding  equation  (F.l)  into  a  Taylor  series  expansion, 
the  resulting  equation  becomes 

p=p  =  l-  ey|-+€2  ^-^(u2  +  *  ^<J>)  -u  +  *5  j 


t  * 


(F.  2) 


Recall  that  the  steady  state  solution  was  represented  in  equation  (2.35)  by 


-\h  e2 


p  =  e 


(6)1. 


When  expanding  (F.3)  into  its  Taylor  series  expansion,  the  result  becomes 


p  =  1  -  he2 


(F.4) 
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where  p  is  the  steady  state  acoustic  pressure.  Therefore ,  the  difference 
in  general  acoustic  pressure  and  steady  state  pressure  can  be  expressed 
by  subtracting  equation  (F.u)  from  (F.2).  For  this  investigation!  a 
restriction  cn  the  velocity  potential  $  was  that  it  was  a  function  of  8 
and  t  only.  In  doing  this(  the  pressure  difference  equation  becomes 


P  -  P 


_e  li  +  siL(l 
e  3t  +  IT  [l? 


2  .P'2' 


(F.5) 


Using  the  same  Fourier  series  expansion  for  the  velocity  potential  ♦  as 
expressed  in  equation  (3.18),  the  acoustic  pressure  difference  equation 
(F.5)  can  be  expressed  in  terms  of  the  product  of  modal  amplitudes  and 
trignometric  function  in  the  transverse  8  direction.  Substituting  the 
appropriate  forms  of  equation  (3.18)  into  equation  (F.5)  and  simplyfying, 
the  resulting  pressure  difference  equation  become 


-  f  df  r  /  dg  dg 

E-l-E.  .  +  e  [-<fxf2  +  glg2)  t  h  ^ 


( dg  dg^  df  df  V 
"  *  df2-  dt7. 


;  Jt-* c  '  ci)  * 


COS0 


cos20 


82  "  f2gl)  +  laldtJ 


fdf  dg  dg  df 


dt  "  JtJ 


sinO 


sin20  . 


(F.6) 
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Since  the  coefficients  in  equation  (F.6)  are  functions  of  time  on.ly, 

\ 

these  coefficients  have  been  included  in  the  calculations  of  the  program 
in  Appendix  B.  Thus,  fcr  any  given  angle  6,  values  for  the  modal 
amplitude  at  any  given  time  range  can  be  calculated  therefore  determining 
the  acoustic  pressure  difference  of  that  desired  location. 


Q 
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